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Abstract 

Randomized algorithms for very large matrix problems have received a great deal of attention in 
recent years. Much of this work was motivated by problems in large-scale data analysis, largely since 
matrices are popular structures with which to model data drawn from a wide range of application 
domains, and this work was performed by individuals from many different research communities. 
While the most obvious benefit of randomization is that it can lead to faster algorithms, either in 
worst-case asymptotic theory and/or numerical implementation, there are numerous other benefits 
that are at least as important. For example, the use of randomization can lead to simpler algorithms 
that are easier to analyze or reason about when applied in counterintuitive settings; it can lead to 
algorithms with more interpretable output, which is of interest in applications where analyst time 
rather than just computational time is of interest; it can lead implicitly to regularization and more 
robust output; and randomized algorithms can often be organized to exploit modern computational 
architectures better than classical numerical methods. 

This monograph will provide a detailed overview of recent work on the theory of randomized matrix 
algorithms as well as the application of those ideas to the solution of practical problems in large-scale 
data analysis. Throughout this review, an emphasis will be placed on a few simple core ideas that 
underlie not only recent theoretical advances but also the usefulness of these tools in large-scale data 
applications. Crucial in this context is the connection with the concept of statistical leverage. This 
concept has long been used in statistical regression diagnostics to identify outliers; and it has recently 
proved crucial in the development of improved worst-case matrix algorithms that are also amenable 
to high-quality numerical implementation and that are useful to domain scientists. This connection 
arises naturally when one explicitly decouples the effect of randomization in these matrix algorithms 
from the underlying linear algebraic structure. This decoupling also permits much finer control in 
the application of randomization, as well as the easier exploitation of domain knowledge. 

Most of the review will focus on random sampling algorithms and random projection algorithms 
for versions of the linear least-squares problem and the low-rank matrix approximation problem. 
These two problems are fundamental in theory and ubiquitous in practice. Randomized methods 
solve these problems by constructing and operating on a randomized sketch of the input matrix A — 
for random sampling methods, the sketch consists of a small number of carefully-sampled and rescaled 
columns/rows of A, while for random projection methods, the sketch consists of a small number of 
linear combinations of the columns/rows of A. Depending on the specifics of the situation, when com- 
pared with the best previously-existing deterministic algorithms, the resulting randomized algorithms 
have worst-case running time that is asymptotically faster; their numerical implementations are faster 
in terms of clock-time; or they can be implemented in parallel computing environments where existing 
numerical algorithms fail to run at all. Numerous examples illustrating these observations will be 
described in detail. 



'Version appearing as a monograph in Now Publishers' "Foundations and Trends in Machine Learning" series. 
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1 Introduction 

This monograph wiU provide a detailed overview of recent work on the theory of randomized matrix 
algorithms as well as the application of those ideas to the solution of practical problems in large-scale 
data analysis. By "randomized matrix algorithms," we refer to a class of recently-developed random 
sampling and random projection algorithms for ubiquitous linear algebra problems such as least-squares 
regression and low-rank matrix approximation. These and related problems are ubiquitous since matrices 
are fundamental mathematical structures for representing data drawn from a wide range of application 
domains. Moreover, the widespread interest in randomized algorithms for these problems arose due to 
the need for principled algorithms to deal with the increasing size and complexity of data that are being 
generated in many of these application areas. 

Not surprisingly, algorithmic procedures for working with matrix-based data have been developed from 
a range of diverse perspectives by researchers from a wide range of areas — including, e.g., researchers from 
theoretical computer science (TCS), numerical linear algebra (NLA), statistics, applied mathematics, data 
analysis, and machine learning, as well as domain scientists in physical and biological sciences — and in 
many of these cases they have drawn strength from their domain-specific insight. Although this has been 
great for the development of the area, and for the "technology transfer" of theoretical ideas to practical 
applications, the technical aspects of dealing with any one of those areas has obscured for many the 
simplicity and generality of some of the underlying ideas; thus leading researchers to fail to appreciate 
the underlying connections and the significance of contributions by researchers outside their own area. 
Thus, rather than focusing on the technical details of proving worst-case bounds or of providing high- 
quality numerical implementations or of relating to traditional machine learning tools or of using these 
algorithms in a particular physical or biological domain, in this review we will focus on highlighting for a 
broad audience the simplicity and generality of some core ideas — ideas that are often obscured but that 
are fruitful for using these randomized algorithms in large-scale data applications. To do so, we will focus 
on two fundamental and ubiquitous matrix problems — least-squares approximation and low-rank matrix 
approximation — that have been at the center of these recent developments. 

The work we will review here had its origins within TCS. In this area, one typically considers a 
particular well-defined problem, and the goal is to prove bounds on the running time and quality-of- 
approximation guarantees for algorithms for that particular problem that hold for "worst-case" input. 
That is, the bounds should hold for any input matrix, independent of any "niceness" assumptions such 
as, e.g., that the elements of the matrix satisfy some smoothness or normalization condition or that 
the spectrum of the matrix satisfies some decay condition. Clearly, the generality of this approach 
means that the bounds will be suboptimal — and thus can be improved — in any particular application 
where stronger assumptions can be made about the input. Importantly, though, it also means that the 
underlying algorithms and techniques will be broadly applicable even in situations where such assumptions 
do not apply. 

An important feature in the use of randomized algorithms in TCS more generally is that one must 
identify and then algorithmically deal with relevant "non-uniformity structure" in the dataQ For the 
randomized matrix algorithms to be reviewed here and that have proven useful recently in NLA and 
large-scale data analysis applications, the relevant non-uniformity structure is defined by the so-called 
statistical leverage scores. Defined more precisely below, these leverage scores are basically the diagonal 
elements of the projection matrix onto the dominant part of the spectrum of the input matrix. As such. 



^For example, for those readers familiar with Markov chain-based Monte Carlo algorithms as used in statistical physics, 
this non-uniformity structure is given by the Boltzmann distribution, in which case the algorithmic question is how to sample 
efficiently with respect to it as an importance sampling distribution without computing the intractable partition function. 
Of course, if the data are sufficiently nice (or if they have been sufficiently preprocessed, or if sufficiently strong assumptions 
are made about them, etc.), then that non-uniformity structure might be uniform, in which case simple methods like uniform 
sampling might be appropriate — but this is far from true in general, either in worst-case theory or in practical applications. 



4 



M. W. Mahoney 



they have a long history in statistical data analysis, where they have been used for outlier detection in 
regression diagnostics. More generally, and very importantly for practical large-scale data applications of 
recently-developed randomized matrix algorithms, these scores often have a very natural interpretation 
in terms of the data and processes generating the data. For example, they can be interpreted in terms 
of the leverage or influence that a given data point has on, say, the best low-rank matrix approximation; 
and this often has an interpretation in terms of high-degree nodes in data graphs, very small clusters in 
noisy data, coherence of information, articulation points between clusters, etc. 

Historically, although the first generation of randomized matrix algorithms (to be described in Sec- 
tion [3]) achieved what is known as additive-error bounds and were extremely fast, requiring just a few 
passes over the data from external storage, these algorithms did not gain a foothold in NLA and only 
heuristic variants of them were used in machine learning and data analysis applications. In order to 
"bridge the gap" between NLA, TCS, and data applications, much finer control over the random sampling 
process was needed. Thus, in the second generation of randomized matrix algorithms (to be described 
in Sections |4] and [5|) that has led to high-quality numerical implementations and useful machine learning 
and data analysis applications, two key developments were crucial. 

• Decoupling the randomization from the linear algebra. This was originally implicit within 
the analysis of the second generation of randomized matrix algorithms, and then it was made 
explicit. By making this decoupling explicit, not only were improved quality-of-approximation 
bounds achieved, but also much finer control was achieved in the application of randomization. For 
example, it permitted easier exploitation of domain expertise, in both numerical analysis and data 
analysis applications. 

• Importance of statistical leverage scores. Although these scores have been used historically 
for outlier detection in statistical regression diagnostics, they have also been crucial in the recent 
development of randomized matrix algorithms. Roughly, the best random sampling algorithms use 
these scores to construct an importance sampling distribution to sample with respect to; and the 
best random projection algorithms rotate to a basis where these scores are approximately uniform 
and thus in which uniform sampling is appropriate. 

As will become clear, these two developments are very related. For example, once the randomization 
was decoupled from the linear algebra, it became nearly obvious that the "right" importance sampling 
probabilities to use in random sampling algorithms are those given by the statistical leverage scores, 
and it became clear how to improve the analysis and numerical implementation of random projection 
algorithms. It is remarkable, though, that statistical leverage scores define the non- uniformity structure 
that is relevant not only to obtain the strongest worst-case bounds, but also to lead to high-quality 
numerical implementations (by numerical analysts) as well as algorithms that are useful in downstream 
scientific applications (by machine learners and data analysts). 

Most of this review will focus on random sampling algorithms and random projection algorithms for 
versions of the linear least-squares problem and the low-rank matrix approximation problem. Here is a 
brief summary of some of the highlights of what follows. 

• Least-squares approximation. Given an m x n matrix A, with m ^ n, and an m-dimensional 
vector 5, the overconstrained least-squares approximation problem looks for the vector Xopt — 
argmin^ll^la; — b\\2. This problem typically arises in statistical models where the rows of A and 
elements of b correspond to constraints and the columns of A and elements of x correspond to vari- 
ables. Classical methods, including the Cholesky decomposition, versions of the QR decomposition, 
and the Singular Value Decomposition, compute a solution in 0{mn^) time. Randomized methods 
solve this problem by constructing a randomized sketch of the matrix A — for random sampling 
methods, the sketch consists of a small number of carefully-sampled and rescaled rows of A (and 
the corresponding elements of b), while for random projection methods, the sketch consists of a 
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small number of linear combinations of the rows of A and elements of b. If one then solves the (still 
overconstrained) subproblem induced on the sketch, then very fine relative-error approximations to 
the solution of the original problem are obtained. In addition, for a wide range of values of m and 
n, the running time is o{'mn^) — for random sampling algorithms, the computational bottleneck is 
computing appropriate importance sampling probabilities, while for random projection algorithms, 
the computational bottleneck is implementing the random projection operation. Alternatively, if 
one uses the sketch to compute a preconditioner for the original problem, then very high-precision 
approximations can be obtained by then calling classical numerical iterative algorithms. Depend- 
ing on the specifics of the situation, these numerical implementations run in o(mn^) time; they are 
faster in terms of clock-time than the best previously-existing deterministic numerical implemen- 
tations; or they can be implemented in parallel computing environments where existing numerical 
algorithms fail to run at all. 

• Low-rank matrix approximation. Given an m x n matrix A and a rank parameter k, the low- 
rank matrix approximation problem is to find a good approximation to A of rank k <C min{m, n}. 
The Singular Value Decomposition provides the best rank-fc approximation to A, in the sense that 
by projecting A onto its top k left or right singular vectors, then one obtains the best approximation 
to A with respect to the spectral and Frobenius norms. The running time for classical low-rank 
matrix approximation algorithms depends strongly on the specifics of the situation — for dense ma- 
trices, the running time is typically O(mnfc); while for sparse matrices, classical Krylov subspace 
methods are used. As with the least-squares problem, randomized methods for the low-rank matrix 
approximation problem construct a randomized sketch — consisting of a small number of either ac- 
tual columns or linear combinations of columns — of the input A, and then this sketch is manipulated 
depending on the specifics of the situation. For example, random sampling methods can use the 
sketch directly to construct relative-error low-rank approximations such as CUR decompositions 
that approximate A based on a small number of actual columns of the input matrix. Alternatively, 
random projection methods can improve the running time for dense problems to O(TOnlogfc); and 
while they only match the running time for classical methods on sparse matrices, they lead to more 
robust algorithms that can be reorganized to exploit parallel computing architectures. 

These two problems are the main focus of this review since they are both fundamental in theory and 
ubiquitous in practice and since in both cases novel theoretical ideas have already yielded practical results. 
Although not the main focus of this review, other related matrix-based problems to which randomized 
methods have been applied will be referenced at appropriate points. 

Clearly, when a very new paradigm is compared with very well-established methods, a nai've imple- 
mentation of the new ideas will perform poorly by traditional metrics. Thus, in both data analysis and 
numerical analysis applications of this randomized matrix algorithm paradigm, the best results have been 
achieved when coupling closely with more traditional methods. For example, in data analysis applica- 
tions, this has meant working closely with geneticists and other domain experts to understand how the 
non-uniformity structure in the data is useful for their downstream applications. Similarly, in scientific 
computation applications, this has meant coupling with traditional numerical methods for improving 
quantities like condition numbers and convergence rates. When coupling in this manner, however, qual- 
itatively improved results have already been achieved. For example, in their empirical evaluation of the 
random projection algorithm for the least-squares approximation problem, to be described in Sections 14.41 
and l4.5l below. Avron, Maymounkov, and Toledo [1] began by observing that "Randomization is arguably 
the most exciting and innovative idea to have hit linear algebra in a long time;" and since their implemen- 
tation "beats LAPACK'fH direct dense least-squares solver by a large margin on essentially any dense tall 



"^Lapack (short for Linear Algebra PACKage) is a high-quality and widely- used software library of numerical routines 
for solving a wide range of numerical linear algebra problems. 



6 



M. W. Mahoney 



matrix," they concluded that their empirical results "show the potential of random sampling algorithms 
and suggest that random projection algorithms should be incorporated into future versions of Lapack." 

The remainder of this review will cover these topics in greater detail. To do so, we will start in 
Section [2] with a few motivating applications from one scientific domain where these randomized matrix 
algorithms have already found application, and we will describe in Section [3] general background on 
randomized matrix algorithms, including precursors to those that are the main subject of this review. 
Then, in the next two sections, we will describe randomized matrix algorithms for two fundamental 
matrix problems: Section |4] will be devoted to describing several related algorithms for the least-squares 
approximation problem; and Section [5] will be devoted to describing several related algorithms for the 
problem of low-rank matrix approximation. Then, Section [6] will describe in more detail some of these 
issues from an empirical perspective, with an emphasis on the ways that statistical leverage scores have 
been used more generally in large-scale data analysis; Section [7] will provide some more general thought 
on this successful technology transfer experience; and Section |8] will provide a brief conclusion. 

2 Matrices in large-scale scientific data analysis 

In this section, we will provide a brief overview of examples of applications of randomized matrix algo- 
rithms in large-scale scientific data analysis. Although far from exhaustive, these examples should serve 
as a motivation to illustrate several complementary perspectives that one can adopt on these techniques. 

2.1 A brief background 

Matrices arise in machine learning and modern massive data set (MMDS) analysis in many guises. One 
broad class of matrices to which randomized algorithms have been applied is object- feature matrices. 

• Matrices from object- feature data. An mxn real- valued matrix A provides a natural structure 
for encoding information about m objects, each of which is described by n features. In astronomy, 
for example, very small angular regions of the sky imaged at a range of electromagnetic frequency 
bands can be represented as a matrix — in that case, an object is a region and the features are the 
elements of the frequency bands. Similarly, in genetics, DNA Single Nucleotide Polymorphism or 
DNA microarray expression data can be represented in such a framework, with Aij representing the 
expression level of the i-th gene or SNP in the j-th experimental condition or individual. Similarly, 
term-document matrices can be constructed in many Internet applications, with Aij indicating the 
frequency of the j-th term in the i-th document. 

Matrices arise in many other contexts — e.g., they arise when solving partial differential equations in 
scientific computation as discretizations of continuum operators; and they arise as so-called kernels when 
describing pairwise relationships between data points in machine learning. In many of these cases, certain 
conditions — e.g., that the spectrum decays fairly quickly or that the matrix is structured such that it 
can be applied quickly to arbitrary vectors or that the elements of the matrix satisfy some smoothness 
conditions — are known or are thought to hold. 

A fundamental property of matrices that is of broad applicability in both data analysis and scientific 
computing is the Singular Value Decomposition (SVD). li A G K™^", then there exist orthogonal matrices 
U = [u^u^ . . . u™] G R'"^^™ and V = [v^v^ ...v"]e E"^" such that U^AV = S = diag(CTi, . . . , a^), where 
E S R™^", v = min{rn, n} and ai > o'2 > • ■ • > > 0. The ai are the singular values of A, the column 
vectors w*, are the i-th left and the i-th right singular vectors of A, respectively. If fc < r = rank(yl), 
then the SVD of A may be written as 

A =. U^V^ - [ Uk Uk,± ] 











(1) 
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Here, is the k x k diagonal matrix containing the top k singular values of A, and 'Sk,± is the {r — k) x 
(r — k) diagonal matrix containing the bottom r — k nonzero singular values of A, is the kx n matrix 
consisting of the corresponding top k right singular vectors!^ etc. By keeping just the top k singular 
vectors, the matrix Ai^ = f/^SfeVJ?" is the best rank-fc approximation to A, when measured with respect 
to the spectral and Frobenius norm. Let = denote the square of the Frobenius 

norm; let \\A\\2 = supj.g]{ri ^.-^q ||^a;||2/||x||2 denote the spectral normQ and, for any matrix A G R™^", 
let ^(i), j S [m] denote the i-th row of A as a row vector, and let A'^^^j G [n] denote the j-th column of 
^ as a column vector. 

Finally, since they will play an important role in later developments, the statistical leverage scores of 
an m X n matrix, with m > n, are defined here. 

Definition 1 Given an arbitrary mx n matrix A, with m > n, let U denote the m x n matrix consisting 
of the n left singular vectors of A, and let U(i) denote the i-th row of the matrix U as a row vector. Then, 
the quantities 

(■i = \\U(i)\\l, forie{l,...,m}, 
are the statistical leverage scores of the rows of A. 

Several things are worth noting about this definition. First, although we have defined these quantities 
in terms of a particular basis, they clearly do not depend on that particular basis, but instead only on 
the space spanned by that basis. To see this, let Pa denote the projection matrix onto the span of the 
columns of A; then, £i = ||f7(i)||| = (UU'^)^^ = (-Pa)^^- That is, the statistical leverage scores of a matrix 
A are equal to the diagonal elements of the projection matrix onto the span of its columns. Second, if 
m > n, then O(mn^) time suffices to compute all the statistical leverage scores exactly: simply perform 
the SVD or compute a QR decomposition of A in order to obtain any orthogonal basis for the range 
of A, and then compute the Euclidean norm of the rows of the resulting matrix. Third, one could also 
define leverage scores for the columns of such a matrix A, but clearly those are all equal to one unless 
m < n or A is rank-deficient. Fourth, and more generally, given a rank parameter k, one can define the 
statistical leverage scores relative to the best rank-k approximation to A to be the m diagonal elements of 
the projection matrix onto the span of the columns of A^, the best rank-fc approximation to A. Finally, 
the coherence 7 of the rows of A is 7 = maxjg{i_ „jj.^i, i.e., it is the largest statistical leverage score 
oiA. 

2.2 Motivating scientific applications 

To illustrate a few examples where randomized matrix algorithms have already been applied in scientific 
data analysis, recall that "the human genome" consists of a sequence of roughly 3 billion base pairs on 
23 pairs of chromosomes, roughly 1.5% of which codes for approximately 20,000 - 25,000 proteins. A 
DNA microarray is a device that can be used to measure simultaneously the genome-wide response of 
the protein product of each of these genes for an individual or group of individuals in numerous different 
environmental conditions or disease states. This very coarse measure can, of course, hide the individual 



^In the text, we will sometimes overload notation and use to refer to any k X n orthonormal matrix spanning the 
space spanned by the top-A: right singular vectors (and similarly for Uj^. and the left singular vectors). The reason is that 
this basis is used only to compute the importance sampling probabilities — since those probabilities are proportional to the 
diagonal elements of the projection matrix onto the span of this basis, the particular basis does not matter. 

* Since the spectral norm is the largest singular value of the matrix, it is an "extremal" norm in that it measures the 
worst-case stretch of the matrix, while the Frobenius norm is more of an "averaging" norm, since it involves a sum over 
every singular direction. The former is of greater interest in scientific computing and NLA, where one is interested in actual 
columns for the subspaces they define and for their good numerical properties, while the latter is of greater interest in data 
analysis and machine learning, where one is more interested in actual columns for the features they define. Both are of 
interest in this review. 
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differences or polymorphic variations. There are numerous types of polymorphic variation, but the most 
amenable to large-scale applications is the analysis of Single Nucleotide Polymorphisms (SNPs), which 
are known locations in the human genome where two alternate nucleotide bases (or alleles, out of A, C, 
G, and T) are observed in a non-negligible fraction of the population. These SNPs occur quite frequently, 
roughly 1 base pair per thousand (depending on the minor allele frequency), and thus they are effective 
genomic markers for the tracking of disease genes {i.e., they can be used to perform classification into sick 
and not sick) as well as population histories (ie, they can be used to infer properties about population 
genetics and human evolutionary history). 

In both cases, m x n matrices A naturally arise, either as a people-by-gene matrix, in which Aij 
encodes information about the response of the j*'' gene in the i*'' individual/condition, or as people- by- 
SNP matrices, in which Aij encodes information about the value of the j*'' SNP in the individual. 
Thus, matrix computations have received attention in these genetics applications [U [3j IH [5l [6l [7] . To 
give a rough sense of the sizes involved, if the matrix is constructed in the nai've way based on data from 
the International HapMap Project [H |^, then it is of size roughly 400 people by 10^ SNPs, although 
more recent technological developments have increased the number of SNPs to well into the millions 
and the number of people into the thousands and tens-of-thousands. Depending on the size of the data 
and the genetics problem under consideration, randomized algorithms can be useful in one or more of 
several ways. 

For example, a common genetics challenge is to determine whether there is any evidence that the 
samples in the data are from a population that is structured, i.e., are the individuals from a homoge- 
neous population or from a population containing genetically distinct subgroups? In medical genetics, 
this arises in case-control studies, where uncorrected population structure can induce false positives; and 
in population genetics, it arises where understanding the structure is important for uncovering the de- 
mographic history of the population under study. To address this question, it is common to perform a 
procedure such as the following. Given an appropriately-normalized (where, of course, the normalization 
depends crucially on domain-specific considerations) m x n matrix A: 

• Compute a full or partial SVD or perform a QR decomposition, thereby computing the eigenvectors 
and eigenvalues of the correlation matrix AA^ . 

• Appeal to a statistical model selection criterior^ to determine either the number k of principal 
components to keep in order to project the data onto or whether to keep an additional principal 
component as significant. 

Although this procedure could be applied to any data set A, to obtain meaningful genetic conclusions 
one must deal with numerous issues^ In any case, however, the computational bottleneck is typically 
computing the SVD or a QR decomposition. For small to medium-sized data, this is not a problem — 
simply call Matlab or call appropriate routines from Lapack directly. The computation of the full 
eigendecomposition takes 0(min{mn^, m^n}) time, and if only k components of the eigendecomposition 
are needed then the running time is typically 0{mnk) time. (This "typically" is awkward from the 
perspective of worst-case analysis, but it is not usually a problem in practice. Of course, one could 
compute the full SVD in 0(min{mn^, m^n}) time and truncate to obtain the partial SVD. Alternatively, 
one could use a Krylov subspace method to compute the partial SVD in 0{mnk) time, but these methods 
can be less robust. Alternatively, one could perform a rank-revealing QR factorization such as that of Gu 
and Eisenstat TS" and then post-process the factors to obtain a partial SVD. The cost of computing the 



^For example, the model selection rule could compare the top part of the spectrum of the data matrix to that of a 
random matrix of the same size |10l 1111 ; or it could use the full spectrum to compute a test statistic to determine whether 
there is more structure in the data matrix than would be present in a random matrix of the same size [3 1121 . 

®For example, how to normalize the data, how to deal with missing data, how to correct for linkage disequilibrium (or 
correlational) structure in the genome, how to correct for closely-related individuals within the sample, etc. 
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QR decomposition is typically 0{mnk) time, although these methods can require slightly longer time in 
rare cases |T3]. See [M] for a discussion of these topics.) 

Thus, these traditional methods can be quite fast even for very large data if one of the dimensions is 
small, e.g., 10^ individuals typed at 10^ SNPs. On the other hand, if both m and n are large, e.g., 10"^ 
individuals at 10^ SNPs, or 10* individuals at 10^ SNPs, then, for interesting values of the rank parameter 
k, the 0{mnk) running time of even the QR decomposition can be prohibitive. As we will see below, 
however, by exploiting randomness inside the algorithm, one can obtain an 0(mn log fc) running time. 
(All of this assumes that the data matrix is dense and fits in memory, as is typical in SNP applications. 
More generally, randomized matrix algorithms to be reviewed below also help in other computational 
environments, e.g., for sparse input matrices, for matrices too large to fit into main memory, when one 
wants to reorganize the steps of the computations to exploit modern multi-processor architectures, etc. 
See [14] for a discussion of these topics.) Since interesting values for k are often in the hundreds, this 
improvement from 0{k) to 0(log k) can be quite significant in practice; and thus one can apply the above 
procedure for identifying structure in DNA SNP data on much larger data sets than would have been 
possible with traditional deterministic methods [15] . 

More generally, a common modus operandi in applying NLA and matrix techniques such as PCA and 
the SVD to DNA microarray, DNA SNPs, and other data problems is: 

• Model the pcople-by-gene or people-by-SNP data as an m x n matrix A. 

• Perform the SVD (or related eigen-methods such as PCA or recently-popular manifold- based meth- 
ods [m [171 [18] that boil down to the SVD, in that they perform the SVD or an eigendecomposition 
on nontrivial matrices constructed from the data) to compute a small number of eigengenes or 
eigenSNPs or eigenpeople that capture most of the information in the data matrix. 

• Interpret the top eigenvectors as meaningful in terms of underlying biological processes; or apply 
a heuristic to obtain actual genes or actual SNPs from the corresponding eigenvectors in order to 
obtain such an interpretation. 

In certain cases, such reification may lead to insight and such heuristics may be justified. For instance, if 
the data happen to be drawn from a Guassian distribution, as in Figure lA, then the eigendirections tend 
to correspond to the axes of the corresponding ellipsoid, and there are many vectors that, up to noise, 
point along those directions. In most cases, however, e.g., when the data are drawn from the union of 
two normals (or mixture of two Gaussians), as in Figure IB, such reification is not valid. In general, the 
justification for interpretation comes from domain knowledge and not the mathematics [191 [3l [20] [21]. The 
reason is that the eigenvectors themselves, being mathematically defined abstractions, can be calculated 
for any data matrix and thus are not easily understandable in terms of processes generating the data: 
eigenSNPs (being linear combinations of SNPs) cannot be assayed; nor can eigengenes (being linear 
combinations of genes) be isolated and purified; nor is one typically interested in how eigenpatients 
(being linear combinations of patients) respond to treatment when one visits a physician. 

For this and other reasons, a common task in genetics and other areas of data analysis is the fol- 
lowing: given an input data matrix A and a parameter k, find the best subset of exactly k actual DNA 
SNPs or actual genes, i.e., actual columns or rows from A, to use to cluster individuals, reconstruct 
biochemical pathways, reconstruct signal, perform classification or inference, etc. Unfortunately, com- 
mon formalizations of this algorithmic problem — including looking for the k actual columns that capture 
the largest amount of information or variance in the data or that are maximally uncorrelated — lead to 
intractable optimization problems [22l[23]. For example, consider the so-called Column Subset Selection 
Problem [24j : given as input an arbitrary m x n matrix A and a rank parameter k, choose the set of 
exactly k columns of A s.t. the m x k matrix C minimizes (over all (^) sets of such columns) the error: 



minWA- PcA\\„ ^inm\\A-CC+A\\ 



(2) 
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where v e {2, F} represents the spectral or Frobenius norra of A, is the Moore-Penrose pseudoinverse 
of C, and Pc — CC'^ is the projection onto the subspace spanned by the columns of C. As we will see 
below, however, by exploiting randomness inside the algorithm, one can find a small set of actual columns 
that is provably nearly optimal. Moreover, this algorithm and obvious heuristics motivated by it have 
already been applied successfully to problems of interest to geneticists such as genotype reconstruction in 
unassayed populations, identifying substructure in heterogeneous populations, and inference of individual 
ancestry [^1 [IS [H [26l [13 [28l [29] . 

In order to understand better the reification issues in scientific data analysis, consider a synthetic data 
set — it was originally introduced in [30] to model oscillatory and exponentially decaying patterns of gene 
expression from [31], although it could just as easily be used to describe oscillatory and exponentially 
decaying patterns in stellar spectra, etc. The data matrix consists of 14 expression level assays (columns 
of A) and 2000 genes (rows of A), corresponding to a 2000 x 14 matrix A. Genes have one of three types of 
transcriptional response: noise (1600 genes); noisy sine pattern (200 genes); and noisy exponential pattern 
(200 genes). Figures IC and ID present the "biological" data, i.e., overlays of five noisy sine wave genes 
and five noisy exponential genes, respectively; Figure IE presents the first and second singular vectors of 
the data matrix, along with the original sine pattern and exponential pattern that generated the data; 
and Figure IF shows that the data cluster well in the space spanned by the top two singular vectors, 
which in this case account for 64% of the variance in the data. Note, though, that the top two singular 
vectors both display a linear combination of oscillatory and decaying properties; and thus they are not 
easily interpretable as "latent factors" or "fundamental modes" of the original (sinusoid and exponential) 
"biological" processes generating the data. This is problematic more generally when one is interested in 
extracting insight or "discovering knowledge" from the output of data analysis algorithms [H] 

Broadly similar issues arise in many other MMDS (modern massive data sets) application areas. In 
astronomy, for example, PCA and the SVD have been used directly for spectral classification [32l [33l 
[34l l35] , to predict morphological types using galaxy spectra [36] , to select quasar candidates from sky 
surveys [37], etc. [38] |39l HO] [41] . Size is an issue, but so too is understanding the data [42] [43]; and many 
of these studies have found that principal components of galaxy spectra (and their elements) correlate 
with various physical processes such as star formation (via absorption and emission line strengths of, 
e.^., the so-called Ha spectral line) as well as with galaxy color and morphology. In addition, there are 
many applications in scientific computing where low-rank matrices appear, e.j;., fast numerical algorithms 
for solving partial differential equations and evaluating potential fields rely on low-rank approximation 
of continuum operators |441 I45| . and techniques for model reduction or coarse graining of multiscale 
physical models that involve rapidly oscillating coeflScients often employ low-rank linear mappings [i5] . 
Recent work that has already used randomized low-rank matrix approximations based on those reviewed 
here include [IT] US] HO] [50] [ST] [5^ • More generally, many of the machine learning and data analysis 
applications cited below use these algorithms and/or greedy or heuristic variants of these algorithms for 
problems in diagnostic data analysis and for unsupervised feature selection for classification and clustering 
problems. 

2.3 Randomization as a resource 

The examples described in the previous subsection illustrate two common reasons for using randomization 
in the design of matrix algorithms for large-scale data problems: 

• Faster Algorithms. In some computation-bound applications, one simply wants faster algorithms 



'^Indeed, after describing the many uses of the vectors provided by the SVD and PCA in DNA microarray analysis, 
Kuruvilla et al. [3] bluntly conclude that "While very efficient basis vectors, the (singular) vectors themselves are completely 
artificial and do not correspond to actual (DNA expression) profiles. ... Thus, it would be interesting to try to find basis 
vectors for all experiment vectors, using actual experiment vectors and not artificial bases that offer little insight." 
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Figure 1: Applying the SVD to data matrices A. (A) 1000 points on the plane, corresponding to a 1000 x 2 
matrix A, (and the two principal components) drawn from a multivariate normal distribution. (B) 1000 
points on the plane (and the two principal components) drawn from a more complex distribution, in this 
case the union of two multivariate normal distributions. (C-F) A synthetic data set considered in [30 to 
model oscillatory and exponentially decaying patterns of gene expression from pH] , as described in the 
text. (C) Overlays of five noisy sine wave genes. (D) Overlays of five noisy exponential genes. (E) The 
first and second singular vectors of the data matrix (which account for 64% of the variance in the data), 
along with the original sine pattern and exponential pattern that generated the data. (F) Projection of 
the synthetic data on its top two singular vectors. Although the data cluster well in the low-dimensional 
space, the top two singular vectors are completely artificial and do not offer insight into the oscillatory 
and exponentially decaying patterns that generated the data. 
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that return more-or-less the exact|f| answer. In many of these applications, one thinks of the rank 
parameter k as the numerical rank of the matrixj^ and thus one wants to choose the error parameter 
e such that the approximation is precise on the order of machine precision. 



• Interpretable Algorithms. In other analyst-bound applications, one wants simpler algorithms 
or more-interpretable output in order to obtain qualitative insight in order to pass to a downstream 
analyst^ In these cases, k is determined according to some domain-determined model selection 
criterion, in which case the difference between ak and ak+i may be small or it may be that ffk+i 3> 
Thus, it is acceptable (or even desirable since there is substantial noise in the data) if e is 
chosen such that the approximation is much less precise. 

Thus, randomization can be viewed as a computational resource to be exploited in order to lead to 
"better" algorithms. Perhaps the most obvious sense of better is faster running time, either in worst- 
case asymptotic theory and/or numerical implementation — we will see below that both of these can be 
achieved. But another sense of better is that the algorithm is more useful or easier to use — e.g., it may 
lead to more interpretable output, which is of interest in many data analysis applications where analyst 
time rather than just computational time is of interest. Of course, there are other senses of better — e.g., 
the use of randomization and approximate computation can lead implicitly to regularization and more 
robust output; randomized algorithms can be organized to exploit modern computational architectures 
better than classical numerical methods; and the use of randomization can lead to simpler algorithms 
that are easier to analyze or reason about when applied in counterintuitive setting fO-but these will not 
be the main focus of this review. 



3 Randomization applied to matrix problems 

Before describing recently-developed randomized algorithms for least-squares approximation and low- 
rank matrix approximation that underlie applications such as those described in Section [21 in this section 
we will provide a brief overview of the immediate precursor^^ of that work. 



®Say, for example, that a numerically-stable answer that is precise to, say, 10 digits of significance is more-or-less exact. 
Exact answers are often impossible to compute numerically, in particular for continuous problems, as anyone who has 
studied numerical analysis knows. Although they will not be the main focus of this review, such issues need to be addressed 
to provide high-quality numerical implementations of the randomized algorithms discussed here. 

®Think of the numerical rank of a matrix as being its "true" rank, up to issues associated with machine precision and 
roundoff error. Depending on the application, it can be defined in one of several related ways. For example, if = min{m, n}, 
then, given a tolerance parameter e, one way to define it is the largest k such that (T^_fe_|_i > e • (T^ |53I . 

^"The tension between providing more interpretable decompositions versus optimizing any single criterion — say, obtaining 
slightly better running time (in scientific computing) or slightly better prediction accuracy (in machine learning) — is well- 
known |21) . It was illustrated most prominently recently by the Netfiix Prize competition — whereas a half dozen or so base 
models captured the main ideas, the winning model was an ensemble of over 700 base models |54| . 

^^Recall that ai is the i*^ singular value of the data matrix. 

^■^Randomization can also be useful in less obvious ways — e.g., to deal with pivot rule issues in communication-constrained 
linear algebra 1551 . or to achieve improved condition number properties in optimization applications [56,. 

^''Although this "first-generation" of randomized matrix algorithms was extremely fast and came with provable quality-of- 
approximation guarantees, most of these algorithms did not gain a foothold in NLA and only heuristic variants of them were 
used in machine learning and data analysis applications. Understanding them, though, was important in the development 
of a "second-generation" of randomized matrix algorithms that were embraced by those communities. For example, while 
in some cases these first-generation algorithms yield to a more sophisticated analysis and thus can be implemented directly, 
more often these first-generation algorithms represent a set of primitives that are more powerful once the randomness is 
decoupled from the linear algebraic structure. 
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3.1 Random sampling and random projections 

Given an m x n matrix yl, it is often of interest to sample randomly a small number of actual columns from 
that matrix^ (To understand why sampling columns (or rows) from a matrix is of interest, recall that 
matrices are "about" their columns and rows j59j — that is, linear combinations are taken with respect to 
them; one all but understands a given matrix if one understands its column space, row space, and null 
spaces; and understanding the subspace structure of a matrix sheds a great deal of light on the linear 
transformation that the matrix represents.) A nai've way to perform this random sampling would be to 
select those columns uniformly at random in i.i.d. trials. A more sophisticated and much more powerful 
way to do this would be to construct an importance sampling distribution and then perform the 

random sample according to it. 

To illustrate this importance sampling approach in a simple setting, consider the problem of approx- 
imating the product of two matrices. Given as input any arbitrary m x n matrix A and any arbitrary 
n X p matrix B: 

• Compute the importance sampling probabilities {pi}f^i, where 

Il^«ll2||%)|l2 



e:up('''ii2||b 



(3) 



Randomly select (and rescale appropriately — if the j*'' column of A is chosen, then scale it by 
l/y^cpj; see ^60^ for details) c columns of A and the corresponding rows of B (again rescaling in the 
same manner), thereby forming m x c and c x p matrices C and R, respectively. 

Two quick points are in order regarding the sampling process in this and other randomized algorithms 
to be described below. First, the sampling here is with replacement. Thus, in particular, if c = n one 
does not necessarily recover the "exact" answer, but of course one should think of this algorithm as 
being most appropriate when c n. Second, if a given column-row pair is sampled, then it must be 
rescaled by a factor depending on the total number of samples to be drawn and the probability that that 
given column- row pair was chosen. The particular form of 1/ cpj ensures that appropriate estimators are 
unbiased; see [B^ for details. 

This algorithm (as well as other algorithms that sample based on the Euclidean norms of the input 
matrices) requires just two passes over the data from external storage. Thus, it can be implemented in 
pass-efficient |60) or streaming |61) models of computation. This algorithm is described in more detail 
in |60) . where it is shown that Frobenius norm bounds of the form 

\\AB-CR\\p<^\\A\\f\\B\\f, (4) 

where 0(1) refers to some constant, hold both in expectation and with high probability. (Issues associated 
with potential failure probabilities, big-0 notation, etc. for this pedagogical example are addressed 
in [BD] — these issues will be addressed in more detail for the algorithms of the subsequent sections.) 
Moreover, if, instead of using importance sampling probabilities of the form ([3|) that depend on both A 
and B, one uses probabilities of the form 

= (5) 



Alternatively, one might be interested in sampling other things like elements 1571 or submatrices [581. Like the algorithms 
described in this section, these other sampling algorithms also achieve additive-error bounds. We will not describe them in 
this review since, although of interest in TCS, they have not (yet?) gained traction in either NLA or in machine learning 
and data analysis applications. 
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that depend on only A (or alternatively ones that depend only on B), then (slightly weaker) bounds of 
the form (|4]) still hold [60]. As we will see, this algorithm (or variants of it, as well as their associated 
bounds) is a primitive that underlies many of the randomized matrix algorithms that have been developed 
in recent years; for very recent examples of this, see [62j l63]. 

To gain insight into "why" this algorithm works, recall that the product AB may be written as 
the outer product or sum of n rank one matrices AB = X]r=i ^^*''-^(*)- When matrix multiplication 
is formulated in this manner, a simple randomized algorithm to approximate the product matrix AB 
suggests itself: randomly sample with replacement from the terms in the summation c times, rescale each 
term appropriately, and output the sum of the scaled terms. If m = p = 1 then A^*\B(^t) G M and it 
is straightforward to show that this sampling procedure produces an unbiased estimator for the sum. 
When the terms in the sum are rank one matrices, similar results hold. In either case, using importance 
sampling probabilities to exploit non- uniformity structure in the data — e.g., to bias the sample toward 
"larger" terms in the sum, as ([3]) does — produces an estimate with much better variance properties. For 
example, importance sampling probabilities of the form ([3|) are optimal with respect to minimizing the 
expectation of \\AB — CR\\f. 

The analysis of the Frobenius norm bound @ is quite simple (60) , using very simple linear algebra and 
only elementary probability, and it can be improved. Most relevant for the randomized matrix algorithms 
of this review is the bound of [64l |65] , where much more sophisticated methods were used to shown that 
if i? = is an n X fc orthogonal matrix Q (i.e., its k columns consist of k orthonormal vectors in R")!^ 
then, under certain assumptions satisfied by orthogonal matrices, spectral norm bounds of the form 

||/ - cc^W^ = \\Q^Q - Q^ss^qW^ < 0{l)^|^^ (6) 

hold both in expectation and with high probability. In this and other cases below, one can represent 
the random sampling operation with a random sampling matrix S — e.g., if the random sampling is 
implemented by choosing c columns, one in each of c i.i.d. trials, then the n x c matrix S has entries 
Sij = 1/ -sjcpi if the i*'' column is picked in the j*^ independent trial, and Sij — otherwise — in which 
case C = AS. 

Alternatively, given an to x n matrix A, one might be interested in performing a random projection 
by post-multiplying A by an n x random projection matrix VL, thereby selecting ^ linear combinations 
of the columns of A. There are several ways to construct such a matrix. 

• Johnson and Lindcnstrauss consider an orthogonal projection onto a random ^-dimensional space |66j 
where £ = O(logTO), and [57] considers a projection onto £ random orthogonal vectors. (In both 
cases, as well as below, the obvious scaling factor of \fnjl is needed.) 

• |68| and [69j choose the entries of f2 as independent, spherically-symmetric random vectors, the 
coordinates of which are i i.i.d. Gaussian Af(0, 1) random variables. 

• [70] chooses the entries of n x £ matrix as { — random variables and also shows that a 
constant factor — up to 2/3 — of the entries of can be set to 0. 

• [Til [72] [73] choose f2 = DHP, where Z? is a n x n diagonal matrix, where each Da is drawn 
independently from {— 1, -f 1} with probability l/2;_ffisannxn normalized Hadamard transform 
matrix, defined below; and P is an n x £ random matrix constructed as follows: Pij = with 
probability \ — q, where q — 0(log^(m)/n); and otherwise either Pij is drawn from a Gaussian 



i^In this case, Q^Q = 4, HQHj = 1, and ||Q|||, = k. Thus, the right hand side of would be 0{l)^/WJc. The tighter 
spectral norm bounds of the form ||6ll on the approximate product of two orthogonal matrices can be used to show that all 
the singular values of are nonzero and thus that rank is not lost — a crucial step in relative-error and high-precision 

randomized matrix algorithms. 
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distribution with an appropriate variance, or is drawn independently from ^ — ^/TJIq, +^/TJIq^ , 
each with probability q/2. 

As with random sampling matrices, post-multiplication by the nx£ random projection matrix 17 amounts 
to operating on the columns — in this case, choosing linear combinations of columns; and thus pre- 
multiplying by ii^ amounts to choosing a small number of linear combinations of rows. Note that, 
aside from the original constructions, these randomized linear mappings are not random projections in 
the usual linear algebraic sense; but instead they satisfy certain approximate metric preserving proper- 
ties satisfied by "true" random projections, and they are useful much more generally in algorithm design. 
Vanilla application of such random projections has been used in data analysis and machine learning 
applications for clustering and classification of data [111 |75l [76l [77l [78l [T^ . 

An important technical point is that the last Hadamard-based construction is of particular importance 
for fast implementations (both in theory and in practice) . Recall that the (non-normalized) nx n matrix 
of the Hadamard transform Hn may be defined recursively as 



Hn/2 Hn/2 
Hn/2 -Hn/2 



with i/o = 



-1 +1 
-1 -1 



in which case the n x n normalized matrix of the Hadamard transform, to be denoted by H hereafter, is 
equal to -^Hn- (For readers not familiar with the Hadamard transform, note that it is an orthogonal 
transformation, and one should think of it as a real-valued version of the complex Fourier transform. 
Also, as defined here, n is a power of 2, but variants of this construction exist for other values of n.) 
Importantly, applying the randomized Hadamard transform, i.e., computing the product xDH for any 
vector X € R" takes 0{n\ogn) time (or even O(nlogr) time if only r elements in the transformed vector 
need to be accessed). Applying such a structured random projection was first proposed in [71[ 172] . it was 
first applied in the context of randomized matrix algorithms in |801 181) , and there has been a great deal of 
research in recent years on variants of this basic structured random projection that are better in theory or 
in practice [TllSIllHllMllHllSilllllHllHIllHilHSlinD]- For example, one could choose n = DHS, where 
iS* is a random sampling matrix, as defined above, that represents the operation of uniformly sampling a 
small number of columns from the randomized Hadamard transform. 

Random projection matrices constructed with any of these methods exhibit a number of similarities, 
and the choice of which is appropriate depends on the application — e.g., a random unitary matrix or a 
matrix with i.i.d. Gaussian entries may be the simplest conceptually or provide the strongest bounds; 
for TCS algorithmic applications, one may prefer a construction with i.i.d. Gaussian, {—1,-1-1}, etc. 
entries, or randomized Hadamard methods that are theoretically efficient to implement; for numerical 
implementations, one may prefer i.i.d. Gaussians if one is working with structured matrices that can 
be applied rapidly to arbitrary vectors and/or if one wants to be very aggressive in minimizing the 
oversampling factor needed by the algorithm, while one may prefer fast-Fourier-based methods that 
are better by constant factors than simple Hadamard-based constructions when working with arbitrary 
dense matrices. 

Intuitively, these random projection algorithms work since, if rj(j) is the j*'' column of n, then An'-J^ 
is a random vector in the range of A. Thus if one generates several such vectors, they will be linearly- 
independent (with very high probability, but perhaps poorly conditioned), and so one might hope to get 
a good approximation to the best rank-fc approximation to A by choosing k or slightly more than k such 
vectors. Somewhat more technically, one can prove that these random projection algorithms work by 
establishing variants of the basic Johnson-Lindenstrauss (JL) lemma, which states: 

• Any set of n points in a high-dimensional Euclidean space can be embedded (via the constructed 
random projection) into an ^-dimensional Euclidean space, where i is logarithmic in n and inde- 
pendent of the ambient dimension, such that all the pairwise distances are preserved to within an 
arbitrarily-small multiplicative (or lie) factor [66l[F7l[68l[69l[70l[7T 1 [72 l [75 ] . 
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This result can then be apphed to (2) vectors associated with the columns of A. The most obvious (but 
not necessarily the best) such set of vectors is the rows of the original matrix A, in which case one shows 
that the random variable — A(,j/)r2||2 equals — A(i/)||2 in expectation (which is usually easy 

to show) and that the variance is sufficiently small (which is usually harder to show). 

By now, the relationship between sampling algorithms and projection algorithms should be clear. 
Random sampling algorithms identify a coordinate-based non-uniformity structure, and they use it to 
construct an importance sampling distribution. For these algorithms, the "bad" case is when that distri- 
bution is extremely nonuniform, i.e., when most of the probability mass is localized on a small number of 
columns. This is the bad case for sampling algorithms in the sense that a naive method like uniform sam- 
pling will perform poorly, while using an importance sampling distribution that provides a bias toward 
these columns will perform much better (at preserving distances, angles, subspaces, and other quantities 
of interest). On the other hand, random projections and randomized Hadamard transforms destroy or 
"wash out" or uniformize that coordinate-based non-uniformity structure by rotating to a basis where the 
importance sampling distribution is very delocalized and thus where uniform sampling is nearly optimal 
(but by satisfying the above JL lemma they too preserve metric properties of interest). For readers more 
familiar with Dirac S functions and sinusoidal functions, recall that a similar situation holds — S functions 
are extremely localized, but when they are multiplied by a Fourier transform, they are converted into 
delocalized sinusoids. As we will see, making such structure explicit has numerous benefits. 

3.2 Randomization for large-scale matrix problems 

Consider the following random projection algorithm that was introduced in the context of understanding 
the success of latent semantic analysis [91 . Given an to x n matrix A and a rank parameter k: 

• Construct an n x £, with £ > a\ogin/e^ for some constant a, random projection matrix ft, as in 
the previous subsection. 

• Return B = Aft. 

This algorithm, which amounts to choosing uniformly a small number £ of columns in a randomly rotated 
basis, was introduced in |91) . where it is proven that 

\\A-PB,MF<\\A-Pu,A\\F + e\\A\\p (7) 

holds with high probability. (Here, i?2fe is the best rank-2fc approximation to the matrix B; Ps^k is the 
projection matrix onto this 2fc-dimensional space; and Pu^. is the projection matrix onto Uk, the top k 
left singular vectors of A.) The analysis of this algorithm boils down to the JL ideas of the previous 
subsection applied to the rows of the input matrix A. That is, the error \\A — Pb^^-^Wf boils down to 
the error incurred by the best rank-2fc approximation plus an additional error term. By applying the 
relative-error JL lemma to the rows of the matrix A, the additional error can be shown to be no greater 
than e \\A\\p. 

Next, consider the following random sampling algorithm that was introduced in the context of clus- 
tering large data sets |92) . Given an m x n matrix A and a rank parameter k: 

• Compute the importance sampling probabilities {pijf^i, where pi = H^'^-'lll/ll^lll^- 

• Randomly select and rescale c = 0[k log fc/e^) columns of A according to these probabilities to form 
the matrix C . 

This algorithm was introduced in [92], although a more complex variant of it appeared in [93] ■ The 
original analysis was extended and simplified in 94" , where it is proven that 

\\A-Pc,A\\2 < \\A~PukA\\, + e\\A\\j, and (8) 
\\A~Pc,A\\p < \\A-PukA\\p + e\\A\\p (9) 
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hold with high probabiUty. (Here, Ck is the best rank-fc approximation to the matrix C, and Pc^. 
is the projection matrix onto this fc-dimensional space.) This additive-error column-based matrix de- 
composition, as well as heuristic variants of it, has been applied in a range of data analysis applica- 
tions [951 EH [Ml EH lag. 

Note that, in a theoretical sense, this and related random sampling algorithms that sample with 
respect to the Euclidean norms of the input columns are particularly appropriate for very large-scale 
settings. The reason is that these algorithms can be implemented efficiently in the pass-efficient or 
streaming models of computation, in which the scarce computational resources are the number of passes 
over the data, the additional RAM space required, and the additional time required. See [60l [94] for 
details about this. 

The analysis of this random sampling algorithm boils down to an approximate matrix multiplication 
result, in a sense that will be constructive to consider in some detail. As an intermediate step in the 
proof of the previous results, that was made explicit in [94], it was shown that 

\\A-Pc,A\\l < \\A~Pu,A\\l + 2\\AA'^ -CC^W^ and 
\\A~Pc,A\\l < \\A~Pu,A\\% + 2Vk\\AA^ ^CC^W^. 

These bounds decouple the linear algebra from the randomization in the following sense: they hold for 
any set of columns, i.e., for any matrix C, and the effect of the randomization enters only through the 
"additional error" term. By using pi = H^^^^Hi/ll^lll^ as the importance sampling probabilities, this 
algorithm is effectively saying that the relevant non-uniformity structure in the data is defined by the 
Euclidean norms of the original matrix. (This may be thought to be an appropriate non-uniformity 
structure to identify since, e.g., it is one that can be identified and extracted in two passes over the data 
from external storage.) In doing so, this algorithm can take advantage of ^ to provide additive-error 
bounds. A similar thing was seen in the analysis of the random projection algorithm — since the JL lemma 
was applied directly to the columns of A, additive-error bounds of the form ([T]) were obtained. 

This is an appropriate point to pause to describe different notions of approximation that a matrix 
algorithm might provide. In the theory of algorithms, bounds of the form provided by ([8]) and ([9]) are 
known as additive- error bounds, the reason being that the "additional" error (above and beyond that 
incurred by the SVD) is an additive factor of the form e times the scale ||v4||^. Bounds of this form are 
very different and in general weaker than when the additional error enters as a multiplicative factor, such 
as when the error bounds are of the form \\A — Pc^A\\ < f{m,n,k,r])\\A — Pjj^AW, where /(•) is some 
function and rj represents other parameters of the problem. Bounds of this type are of greatest interest 
when /(•) does not depend on m or n, in which case they are known as a constant-factor bounds, or when 
they depend on m and n only weakly. The strongest bounds are when / = 1 -I- e, for an error parameter 
e, i.e., when the bounds are of the form ||^ — Pc^^H < (1 + e)||yl — P(7^yl||. These relative- error bounds 
are the gold standard, and they provide a much stronger notion of approximation than additive-error or 
weaker multiplicative-error bounds. We will see bounds of all of these forms below. 

One application of these random sampling ideas that deserves special mention is when the input matrix 
A is symmetric positive semi-definite. Such matrices are common in kernel-based machine learning, and 
sampling columns in this context often goes by the name the Nystrom method. Originating in integral 
equation theory, the Nystrom method was introduced into machine learning in [99] and it was analyzed and 
discussed in detail in [lOOj . Applications to large-scale machine learning problems include [lOH 11021 1103j 
and [104|, I105| I106j , and applications in statistics and signal processing include [107[ I108[ I109| I110| lllll 
11121 1113) . As an example, the Nystrom method can be used to provide an approximation to a matrix 
without even looking at the entire matrix — under assumptions on the input matrix, of course, such as 
that the leverage scores are approximately uniform. 
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3.3 A retrospective and a prospective 

Much of the early work in TCS focused on randomly sampling columns according to an importance 
sampling distribution that depended on the Euclidean norm of those columns [93l [92l |60l [Ml I114| l65] . 
This had the advantage of being "fast," in the sense that it could be performed in a small number 
of "passes" over that data from external storage, and also that additive-error quality-of-approximation 
bounds could be proved. On the other hand, this had the disadvantage of being less immediately- 
applicable to scientific computing and large-scale data analysis applications. At root, the reason is that 
these algorithms didn't highlight "interesting" or "relevant" non-uniformity structure, which then led to 
bounds that were rather coarse. For example, columns are easy to normalize and are often normalized 
during data preprocessing. Even when not normalized, column norms can still be uninformative, as in 
heavy-tailed graph datal]f| where they often correlate strongly with simpler statistics such as node degree. 

Relatedly, bounds of the form (jj]) do not exploit the underlying vector space structure. This is 
analogous to how the JL lemma was applied in the analysis of the random projection algorithm — by 
applying the JL lemma to the actual rows of A, as opposed to some other more refined vectors associated 
with the rows of A, the underlying vector space structure was ignored and only coarse additive-error 
bounds were obtained. To obtain improvements and to bridge the gap between TCS, NLA, and data 
applications, much finer bounds that take into account the vector space structure in a more refined way 
were needed. To do so, it helped to identify more refined structural properties that decoupled the random 
matrix ideas from the underlying linear algebraic ideas — understanding this will be central to the next 
two sections. 

Although these newer algorithms identified more refined structural properties, they have the same 
general structure as the original randomized matrix algorithms. Recall that the general structure of the 
algorithms just reviewed is the following. 

• Preprocess the input by: defining a non-uniformity structure over the columns of the input matrix; 
or performing a random projection/rotation to uniformize that structure. 

• Draw a random sample of columns from the input matrix, either using the non-uniformity structure 
as an importance sampling distribution to select actual columns, or selecting columns uniformly at 
random in the rotated basis. 

• Postprocess the sample with a traditional deterministic NLA method. 

In the above algorithms, the preprocessing was very fast, in that the importance sampling distribution 
could be computed by simply passing over the data a few times from external storage; and the post- 
processing consists of just computing the best rank-fc or best rank-2fc approximation to the sample. As 
will become clear below, by making one or both of these steps more sophisticated, very substantially 
improved results can be obtained, both in theory and in practice. This can be accomplished, e.g., by 
using more sophisticated sampling probabilities or coupling the randomness in more sophisticated ways 
with traditional NLA methods, which in some cases will require additional computation. 

4 Randomized algorithms for least-squares approximation 

In this section and the next, we will describe randomized matrix algorithms for the least-squares approxi- 
mation and low-rank approximation problems. The analysis of low-rank matrix approximation algorithms 
described in Section [5] boils down to a randomized approximation algorithm for the least-squares approx- 
imation problem jll5[ [211 [24l 1116) . For this reason, for pedagogical reasons, and due to the fundamental 



^ By heavy-tailed graph, consider a graph — or equivalently the adjacency matrix of such a graph — in which quantities 
such as the degree distribution or eigenvalue distribution decay in a heavy-tailed or power law manner. 
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importance of the least-squares problem more generally, randomized algorithms for the least-squares 
problem will be the topic of this section. 

4.1 Different perspectives on least-squares approximation 

Consider the problem of finding a vector x such that Ax w 6, where the rows of A and elements of b 
correspond to constraints and the columns of A and elements of x correspond to variables. In the very 
overconstrained least-squares approximation problem, where the m x n matrix A has m 3> n^]\ there is 
in general no vector x such that Ax — b, and it is common to quantify "best" by looking for a vector 
Xopt such that the Euclidean norm of the residual error is small, i.e., to solve the least-squares (LS) 
approximation problem 

Xopt = argmin^llAx - 6II2. (10) 

This problem is ubiquitous in applications, where it often arises from fitting the parameters of a model 
to experimental data, and it is central to theory. Moreover, it has a natural statistical interpretation 
as providing the best estimator within a natural class of estimators; and it has a natural geometric 
interpretation as fitting the part of the vector b that resides in the column space of A. From the viewpoint 
of low-rank matrix approximation, this LS problem arises since measuring the error with a Frobenius or 
spectral norm, as in ([2]), amounts to choosing columns that are "good" in a least squares sense 

There are a number of different perspectives one can adopt on this LS problem. Two major perspec- 
tives of central importance in this review are the following. 

• Algorithmic perspective. From an algorithmic perspective, the relevant question is: how long 
does it take to compute a^opt? The answer to this question is that is takes 0{mn^) time |53] . 
This can be accomplished with one of several algorithms — with the Cholesky decomposition (which 
is good if A has full column rank and is very well-conditioned); or with a variant of the QR 
decomposition (which is somewhat slower, but more numerically stable); or by computing the full 
SVD A — UYiV^ (which is often, but certainly not always, overkill, but which can be easier to 
explair0), and letting Xopt = VYi^U'^b. Although these methods differ a great deal in practice and 
in terms of numerical implementation, asymptotically each of these methods takes a constant times 
mn^ time to compute a vector Xopt- Thus, from an algorithmic perspective, a natural next question 
to ask is: can the general LS problem be solved, either exactly or approximately, in o{mv?) time@ 
with no assumptions at all on the input data? 

• Statistical perspective. From a statistical perspective, the relevant question is: when is com- 
puting the Xopt the right thing to do? The answer to this question is that this LS optimization 
is the right problem to solve when the relationship between the "outcomes" and "predictors" is 



^^In this section only, we will assume that m ^ n. 

^^Intuitively, these low-rank approximation algorithms find columns that provide a space that is good in a least-squares 
sense, when compared to the best rank-fc space, at reconstructing every row of the input matrix. Thus, the reason for the 
connection is that the merit function that describes the quality of those algorithms is typically reconstruction error with 
respect to the spectral or Frobenius norm. 

^^The SVD has been described as the "Swiss Army Knife" of NLA |117l . That is, given it, one can do nearly anything one 
wants, but it is almost always overkill, as one rarely if ever needs its full functionality. Nevertheless, for pedagogical reasons, 
since other decompositions are typically better in terms of running time by only constant factors, and since numerically- 
stable algorithms for these latter decompositions can be quite complex, it is convenient to formulate results in terms of the 
SVD and the best rank-fc approximation to the SVD. 

^''Formally, /(n) = o{g(ri)) as n — >■ 00 means that for every positive constant e there exists a constant A'' such that 
foi' ^11 n > N. Informally, it means that f{n) grows more slowly than g{n). Thus, if the running time of 
an algorithm is o(mn^) time, then it is asymptotically faster than any (arbitrarily small) constant times mn^. 
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roughly linear and when the error processes generating the data are "nice" (in the sense that they 
have mean zero, constant variance, are uncorrelated, and are normally distributed; or when we have 
adequate sample size to rely on large sample theory) |118) . Thus, from a statistical perspective, a 
natural next question to ask is: what should one do when the assumptions underlying the use of 
LS methods are not satisfied or are only imperfectly satisfied? 

Of course, there are also other perspectives that one can adopt. For example, from a numerical perspec- 
tive, whether the algorithm is numerically stable, issues of forward versus backward stability, condition 
number issues, and whether the algorithm takes time that is a large or small constant multiplied by 
mui{rniT? ,rn^n} are of paramount importance. 

When adopting the statistical perspective, it is common to check the extent to which the assumptions 
underlying the use of LS have been satisfied. To do so, it is common to assume that b = Ax + e, 
where b is the response, the columns A'*^ are the carriers, and e is a "nice" error processj^ Then 
Xopt = {A'^A)~^A'^b, and thus b — Hb, where the projection matrix onto the column space of A, 

H = A{A^A)-^A^, 

is the so-called hat matrix. It is known that Hij measures the influence or statistical leverage exerted on 
the prediction bi by the observation bj |119[ I118[ I120[ I121[ I122j . Relatedly, if the z*'' diagonal element 
of H is particularly large then the z*'' data point is particularly sensitive or influential in determining 
the best LS fit, thus justifying the interpretation of the elements Ha as statistical leverage scores |21| . 
These leverage scores have been used extensively in classical regression diagnostics to identify potential 
outliers by, e.g., flagging data points with leverage score greater than 2 or 3 times the average value in 
order to be investigated as errors or potential outliers |118j . Moreover, in the context of recent graph 
theory applications, this concept has proven useful under the name of graph resistance |123j : and, for the 
matrix problems considered here, some researchers have used the term coherence to measure the degree 
of non-uniformity of these statistical leverage scores |124[ I125| I126j . 

In order to compute these quantities exactly, recall that if U is any orthogonal matrix spanning the 
column space of A, then H = Pu = UU^ and thus 

Ha — I |C^(i) 1 121 

i.e., the statistical leverage scores equal the Euclidean norm of the rows of any such matrix U |115[ [21] . 
Recall Definition [T] from Section 12.11 (Clearly, the columns of such a matrix U are orthonormal, but 
the rows of U in general are not — they can be uniform if, e.g., U consists of columns from a truncated 
Hadamard matrix; or extremely nonuniform if, e.g., the columns of U come from a truncated identity 
matrix; or anything in between.) More generally, and of interest for the low-rank matrix approximation 
algorithms in Section [SJ the statistical leverage scores relative to the best rank-k approximation to A are 
the diagonal elements of the projection matrix onto the best rank-fc approximation to A. Thus, they can 
be computed from 

{PUk)i'i. = l|C^fc,(j)||2: 

where C/fc,(i) is the i^^ row of any matrix spanning the space spanned by the top k left singular vectors 
of A (and similarly for the right singular subspace if columns rather than rows are of interest). 

In many diagnostic applications, e.g., when one is interested in exploring and understanding the 
data to determine what would be the appropriate computations to perform, the time to compute or 
approximate {Pu)ii or {Puk)ii is not the bottleneck. On the other hand, in cases where this time is the 



^^This is typically done by assuming that the error process e consists of i.i.d. Gaussian entries. As with the construction 
of random projections in Section l3.ll numerous other constructions are possible and will yield similar results. Basically, it 
is important that no one or small number of data points has a particularly large influence on the LS fit, in order to ensure 
that techniques from large-sample theory like measure concentration apply. 
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bottleneck, an algorithm we will describe in Section 14.4.21 will provide very fine approximations to all 
these leverage scores in time that is qualitatively faster than that required to compute them exactly. 

4.2 A simple algorithm for approximating least-squares approximation 

Returning to the algorithmic perspective, consider the following random sampling algorithm for the LS 
approximation problem [1271 1115] . Given a very overconstrained LS problem, where the input matrix A 
and vector b are arbitrary, but m ^ n: 

• Compute the normalized statistical leverage scores {pi}™ j^, i.e., compute pi — ||C^(i)||2/^i where U 
is the m X n matrix consisting of the left singular vectors of Ar-^l 

• Randomly sample and rescal r — 0{n\ogn/e^) constraints, i.e., rows oi A and the corresponding 
elements of b, using these scores as an importance sampling distribution. 

• Solve the induced subproblem Xopt = argmin^. 1 1 5*^2; — S'6||2, where the r x m matrix S represents 
the sampling-and-rescaling operation. 

The induced subproblem can be solved using any appropriate direct or iterative LS solver as a black box. 
For example, one could compute the solution directly via the generalized inverse as Xopt = (SA)'' Sb, 
which would take 0{rn'^) time [53]. Alternatively, one could use iterative methods such as the Conjugate 
Gradient Normal Residual method, which can produce an e-approximation to the optimal solution of 
the sampled problem in 0{k(S A)rn\og{l/ e)) time, where k{SA) is the condition number of SA [53 . As 
stated, this algorithm will compute all the statistical leverage scores exactly, and thus it will not be faster 
than the traditional algorithm for the LS problem — importantly, we will see below how to get around 
this problem. 

Since this overconstrainccF^ LS algorithm samples constraints and not variables, the dimensionality 
of the vector Xopt that solves the subproblem is the same as that of the vector Xopt that solves the original 
problem. The algorithm just presented is described in more detail in |127| 11151 E], where it is shown 
that relative-error bounds of the form 



hold, where k{A) is the condition number of A and where 7 = ||?7?7'^&||2/||6||2 is a parameter defining 
the amount of the mass of b inside the column space of Of course, there is randomization inside 



Stating this in terms of the singular vectors is a convenience, but it can create confusion. In particular, although 
computing the SVD is sufficient, it is by no means necessary — here, U can be any orthogonal matrix spanning the column 
space of A | 21l . Moreover, these probabilities are robust, in that any probabilities that are close to the leverage scores will 
suffice; see 1601 for a discussion of approximately-optimal sampling probabilities. Finally, as we will describe below, these 
probabilities can be approximated quickly, i.e., more rapidly than the time needed to compute a basis exactly, or the matrix 
can be preprocessed quickly to make them nearly uniform. 

^''Recall from the discussion in Section l3.ll that each sampled row should be rescaled by a factor of 1/rpi. Thus, it is 
these sampled-and-rescaled rows that enter into the subproblem that this algorithm constructs and solves. 

^*Not surprisingly, similar ideas apply to underconstrained LS problems, where m <^ n, and where the goal is to compute 
the minimum-length solution. In this case, one randomly samples columns, and one uses a somewhat more complicated 
procedure to construct an approximate solution, for which relative-error bounds also hold. In fact, as we will describe in 
Section 14.4.21 the quantities {pi}^^ can be approximated in o(mn?) time by relating them to an underconstrained LS 
approximation problem and running such a procedure. 

^^We should reemphasize that such relative-error bounds, either on the optimum value of the objective function, as in JTTJ 
and as is more typical in TCS, or on the vector or "certificate" achieving that optimum, as in II12I I and as is of greater 
interest in NLA, provide an extremely strong notion of approximation. 



\\b-Axopt\\2 < {I + £)\\b - AxoptlU and 



(11) 
(12) 
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this algorithm, and it is possible to flip a fair coin "heads" 100 times in a row. Thus, as stated, with 
r = 0(nlogn/e^), the randomized least-squares algorithm just described might fail with a probability 5 
that is no greater than a constant (say 1/2 or 1/10 or 1/100, depending on the (modest) constant hidden 
in the O(-) notation) that is independent of m and n. Of course, using standard methods jl28] . this 
failure probability can easily be improved to be an arbitrarily small 8 failure probability. For example, 
this holds if r = 0(nlog(n) log(l/(5)/e^) in the above algorithm; alternatively, it holds if one repeats the 
above algorithm 0(log(l/(5)) times and keeps the best of the results. 

As an aside, it is one thing for TCS researchers, well-known to be cavalier with respect to constants 
and even polynomial factors, to make such observations about big-0 notation and the failure probability, 
but by now these facts are acknowledged more generally. For example, a recent review of coupling 
randomized low-rank matrix approximation algorithms with traditional NLA methods |14| starts with 
the following observation. "Our experience suggests that many practitioners of scientific computing view 
randomized algorithms as a desperate and final resort. Let us address this concern immediately. Classical 
Monte Carlo methods are highly sensitive to the random number generator and typically produce output 
with low and uncertain accuracy. In contrast, the algorithms discussed herein are relatively insensitive 
to the quality of randomness and produce highly accurate results. The probability of failure is a user- 
specified parameter that can be rendered negligible (say, less than 10^^^) with a nominal impact on the 
computational resources required." 

Finally, it should be emphasized that modulo this failure probability 5 that can be made arbitrarily 
small without adverse effect and an error e that can also be made arbitrarily small, the above algorithm 
(as well as the basic low-rank matrix approximation algorithms of Section [5] that boil down to this 
randomized approximate LS algorithm) returns an answer Xopt that satisfies bounds of the form (jlip 
and (jl2p . independent of any assumptions at all on the input matrices A and b. 

4.3 A basic structural result 

What the above random sampling algorithm highlights is that the "relevant non-uniformity structure" 
that needs to be dealt with in order to solve the LS approximation problem is defined by the statistical 
leverage scores. (Moreover, since the randomized algorithms for low-rank matrix approximation to be 
described in Section [5] boil down to a least-squares problem, an analogous statement is true for these 
low-rank approximation algorithms.) To see "why" this works, it is helpful to identify a deterministic 
structural condition sufficient for relative-error approximation — doing so decouples the linear algebraic 
part of the analysis from the randomized matrix part. This condition, implicit in the analysis of |127[I115] . 
was made explicit in [81j . 

Consider preconditioning or premultiplying the input matrix A and the target vector b with some 
arbitrary matrix Z, and consider the solution to the LS approximation problem 

Xopt = a3:gimii^\\Z{Ax ~b)\\2. (13) 

Thus, for the random sampling algorithm described in Section [4.21 the matrix Z is a carefully-constructed 
data-dependent random sampling matrix, and for the random projection algorithm below it is a data- 
independent random projection, but more generally it could be any arbitrary matrix Z . Recall that the 
SVD oi A\s A = Ua'^aVJ; and, for notational simplicity, let b-^ = U^U^ b denote the part of the right 
hand side vector b lying outside of the column space of A. Then, the following structural condition holds. 

• Structural condition underlying the randomized least-squares algorithm. Under the 
assumption that Z satisfies the following two conditions: 

aL„ (ZUa) > 1/V2; and (14) 
\\UXZ^Zb^\\l<^\\Axopt-b\\l (15) 
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for some e G (0, 1), the solution vector Xopt to the LS approximation problem (jl3p satisfies relative- 
error bounds of the form pA|) and ([T2|) . 

In this condition, the particular constants l/-\/2 and 1/2 clearly don't matter — they have been chosen 
for ease of comparison with [8T]. Also, recall that ||6^||2 = ||^a;opt — h\\2- Several things should be noted 
about these two structural conditions: 

• First, since ai{UA) — 1, for all i £ [n], Condition indicates that the rank of ZUa is the same 
as that of Ua- Note that although Condition (IT4l) only states that af{ZUA) > l/\/2, for all i £ [n], 
for the randomized algorithm of Section l4?2l it will follow that |l — (7f{ZUA}\ < 1 — 2^^/'^, for all 
i £ [n]. Thus, one should think of Condition (fT4| as stating that ZUa is an approximate isometry. 
Thus, this expression can be bounded with the approximate matrix multiplication spectral norm 
bound of 

• Second, since before preprocessing by Z, b-^ = U^Uj^^b is clearly orthogonal to Ua, Condition ([15]) 
simply states that after preprocessing Zh^ remains approximately orthogonal to ZUa- Although 
Condition ([T5|) depends on the right hand side vector b, the randomized algorithm of Section 14.21 
satisfies it without using any information from b. The reason for not needing information from 
h is that the left hand side of Condition (fTSj) is of the form of an approximate product of two 
different matrices — where for the randomized algorithm of Section 14.21 the importance sampling 
probabilities depend only on one of the two matrices — and thus one can apply an approximate 
matrix multiplication bound of the form Q. 

• Third, as the previous two points indicate. Condition ([T4l) and ([TS]) both boil down to the problem 
of approximating the product of two matrices, and thus the algorithmic primitives on approximate 
matrix multiplication from Section [3.11 will be useful, either explicitly or within the analysis. 

It should be emphasized that there is no randomization in these two structural conditions — they are 
deterministic statements about an arbitrary matrix Z that represent a structural condition sufficient 
for relative-error approximation. Of course, if Z happens to be a random matrix, e.g., representing a 
random projection or a random sampling process, then Conditions (fT4|) or (|15p may fail to be satisfied — 
but, conditioned on their being satisfied, the relative-error bounds of the form pT|) and (fT^ follow. Thus, 
the effect of randomization enters only via Z, and it is decoupled from the linear algebraic structure. 

4.4 Making this algorithm fast — in theory 

Given this structural insight, what one does with it depends on the application. In some cases, as 
in certain genetics applications [101 ESI [21] or when solving the Column Subset Selection Problem as 
described in Section 15.21 the time for the computation of the leverage scores is not the bottleneck, and 
thus they may be computed with the traditional procedure. In other cases, as when dealing with not- 
extremely-sparse random matrices or other situations where it is expected or hoped that no single data 
point is particularly influential, it is assumed that the scores are exactly or approximately uniform, in 
which case uniform sampling is appropriate. In general, however, the simple random sampling algorithm 
of Section l4?2l requires the computation of the normalized statistical leverage scores, and thus it runs in 
O^rnn^) time. 

In this subsection, we will describe two ways to speed up the random sampling algorithm of Section|4?2] 
so that it runs in o{mn^) time for arbitrary input. The first way involves preprocessing the input with 
a randomized Hadamard transform and then calling the algorithm of Section 14.21 with uniform sampling 
probabilities. The second way involves computing a quick approximation to the statistical leverage 
scores and then using those approximations as the importance sampling probabilities in the algorithm of 
Section [321 Both of these methods provide fast and accurate algorithms in theory, and straightforward 
extensions of them provide very fast and very accurate algorithms in practice. 
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4.4.1 A fast random projection algorithm for the LS problem 

Consider the following structured random projection algorithm for approximating the solution to the LS 
approximation problem. 

• Prcmultiply A and b with an n x n randomized Hadamard transform HD. 



ing elements from HDb. 

• Solve the induced subproblem Xopt = aTgimii^\\SHDAx — SHDb\\2, where the r x m matrix S 
represents the sampling operation. 

This algorithm, which first preprocesses the input with a structured random projection and then solves 
the induced subproblem, as well as a variant of it that uses the original "fast" Johnson-Lindenstrauss 
transform [7TJ [72l [73] , was presented in [80l [81] , (where a precise statement of r is given and) where it is 
shown that relative-error bounds of the form (fTTj) and ()12|) hold. 

To understand this result, recall that premultiplication by a randomized Hadamard transform is a 
unitary operation and thus does not change the solution; and that from the SVD of A and of HDA it 
follows that Uhda = HDUa- Thus, the "right" importance sampling distribution for the preprocessed 
problem is defined by the diagonal elements of the projection matrix onto the span of HDUa- Importantly, 
application of such a Hadamard transform tends to "uniformize" the leverage scores0 in the sense that 
all the leverage scores associated with Uhda are (up to logarithmic fluctuations) uniform [71] [ST]. Thus, 
uniform sampling probabilities are optimal, up to a logarithmic factor which can be absorbed into the 
sampling complexity. Overall, this relative-error approximation algorithm for the LS problem run in 



when m ^ n. Although the ideas underlying the Fast Fourier Transform have been around and used in 
many applications for a long time |129[ 1130] , they were first applied in the context of randomized matrix 
algorithms only recently [lU [80l [81] . 

The o(m7i^) running time is most interesting when the input matrix to the overconstrained LS problem 
is dense; if the input matrix is sparse, then a more appropriate comparison might have to do with the 
number of nonzero entries. In general, however, random Gaussian projections and randomized Hadamard- 
based transforms tend not to respect sparsity. In some applications, e.g., the algorithm of [131 that 
is described in Section 14.51 and that automatically speeds up on sparse matrices and with fast linear 
operators, as well as several of the randomized algorithms for low-rank approximation to be described 
in Section 15.31 this can be worked around. In general, however, the question of using sparse random 
projections and sparsity-preserving random projections is an active topic of research (87> ,88i l89l 1132] . 

4.4.2 A fast random sampling algorithm for the LS problem 

Next, consider the following algorithm which takes as input an arbitrary m x n matrix A, with m 3> n, 
and which returns as output approximations to all m of the statistical leverage scores of A. 

• Premultiply ^ by a structured random projection, e.g., fii = SHD from Section [3.11 which repre- 
sents uniformly sampling roughly ri = 0(m\ogn/e) rows from a randomized Hadamard transform. 

• Compute the mxr2 matrix X = A{VliA)^Q,2, where 0,2 is an ri x r2 unstructured random projection 
matrix and where the dagger represents the Moore-Penrose pseudoinverse. 



^^As noted above, this is for very much the same reason that a Fourier matrix delocalizcs a localized 5-function; and it 
also holds for the application of an unstructured random orthogonal matrix or random projection. 




rows from HDA and the correspond- 



o{mn ) time [SOI [81]— essentially 




much less than 0(rnn?) 
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• For each i = 1, . . . , m, compute and return — 



This algorithm was introduced in [133 , based on ideas in jl34j 



In |133) . it is proven that 



- h\ < eii for all i^l 



. . , m 



where ii are the statistical leverage scores of Definition [T] That is, this algorithm returns a relative-error 
approximation to every one of the m statistical leverage scores. Moreover, in |133| it is also proven that 
this algorithm runs in o(mn^) time — due to the structured random projection in the first step, the running 
time is basically the same time as that of the fast random projection algorithm described in the previous 
subsectionl^ In addition, given an arbitrary rank parameter k and an arbitrary-sized mxn matrix A, this 
algorithm can be extended to approximate the leverage scores relative to the best rank-Zc approximation 
to A in roughly 0{mnk) time. See |133j for a discussion of the technical issues associated with this. In 
particular, note that the problem of asking for approximations to the leverage scores relative to the best 
rank-fc space of a matrix is ill-posed; and thus one must replace it by computing approximations to the 
leverage scores for some space that is a good approximation to the best rank-fc space. 

Within the analysis, this algorithm for computing rapid approximations to the statistical leverage 
scores of an arbitrary matrix basically boils down to an wnderconstrained LS problem, in which a struc- 
tured random projection is carefully applied, in a manner somewhat analogous to the fast oi^erconstrained 
LS random projection algorithm in the previous subsection. In particular, let A be an m x n matrix, with 
m ^ n, and consider the problem of finding the minimum-length solution to Xopt = argmin^j — &| I2 = 
A'^b. Sampling variables or columns from A can be represented by postmultiplying A hy a n x c 
(with c > m) column-sampling matrix S to construct the (still underconstrained) least-squares prob- 
lem: Xopt = argmin^||A5'S'-^a; — b\\2 = A'^ (AS)'^^ {AS)^b. The second equality follows by inserting 
P4T = A^A^+ to obtain ASS'^ A^ A^+x - b inside the || • II2 and recahing that A+ = A^ A^+ A+ for the 
Moore-Penrose pseudoinverse. If one constructs S by randomly sampling c = 0{{n/e^) log(n/e)) columns 
according to "column- leverage-score" probabilities, i.e., exact or approximate diagonal elements of the 
projection matrix onto the row space of A, then it can be proven that \\xopt — XoptlU ^ e||a;opt||2 holds, 
with high probability. Alternatively, this underconstrained LS problem problem can also be solved with a 
random projection. By using ideas from the previous subsection, one can show that if S instead represents 
a random projection matrix, then by projecting to a low-dimensional space (which takes o{m'^n) time 
with a structured random projection matrix), then relative-error approximation guarantees also hold. 

Thus, one can run the following algorithm for approximating the solution to the general overcon- 
strained LS approximation problem. 

• Run the algorithm of this subsection to obtain numbers for each i = 1, . . . , m, rescaling them to 
form a probability distribution over {!,..., m}. 

• Call the randomized LS algorithm that is described in Section l4?2l except using these numbers ii 
(rather than the exact statistical leverage scores ii) to construct the importance sampling distribu- 
tion over the rows of A. 

That is: approximate the normalized statistical leverage scores in oimv?) time; and then call the random 
sampling algorithm of Section |421 using those approximate scores as the importance sampling distribution. 
Clearly, this combined algorithm provides relative-error guarantees of the form (jlip and (|12|) , and overall 
it runs in o(rarP') time. 



^^Recall that since the coherence of a matrix is equal to the largest leverage score, this algorithm also computes a relative- 
error approximation to the coherence of the matrix in o(irm?) time — which is qualitatively faster than the time needed to 
compute an orthogonal basis spanning the original matrix. 
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4.4.3 Some additional thoughts 

The previous two subsections have shown that the random samphng algorithm of Section 14.21 which 
naively needs O(mn^) to compute the importance sampling distribution, can be sped up in two different 
ways. One can spend o{mii?) to uniformize the sampling probabilities and then sample uniformly; or 
one can spend olrnn?) time to compute approximately the sampling probabilities and then use those 
approximations as an importance sampling distribution. Both procedures take basically the same time, 
which should not be surprising since the approximation of the statistical leverage scores boils down to an 
underconstrained LS problem; and which procedure is preferable in any given situation depends on the 
downstream application. 

Finally, to understand better the relationship between random sampling and random projection algo- 
rithms for the least squares problem, consider the following vanilla random projection algorithm. Given 
a matrix A and vector 6, representing a very overconstrained LS problem, one can: 

• Construct an 0(nlogn/e^) x m random projection matrix i7, where the entries consist (up to 
scaling) of i.i.d. Gaussians or {—1,+!}. 

• Solve the induced subproblem Xopt = argmin^ 1 1 17 Ax — r2fe||2- 

It is relatively-easy to show that this algorithm also satisfies relative-error bounds of the form pT|) 
and mi). Importantly, though, this random projection algorithm requires O(mn^) time to implement. 
For the random sampling algorithm of Section|421 the 0{mn^) computational bottleneck is computing the 
importance sampling probabilities; while for this random projection algorithm the 0{mn^) computational 
bottleneck is actually performing the matrix-matrix multiplication representing the random projection. 
Thus, one can view the o{mn?) random projection algorithm that uses a small number of rows from a 
randomized Hadamard transform in one of two ways: either as a structured approximation to a usual 
random projection matrix with i.i.d. Gaussian or {—1, -1-1} entries that can be rapidly applied to arbitrary 
vectors; or as preprocessing the input with an orthogonal transformation to uniformize the leverage scores 
so that uniform sampling is appropriate. 

4.5 Making this algorithm fast — in practice 

Several "rubber-hits-the-road" issues need to be dealt with in order for the algorithms of the previous 
subsection to yield to high-precision numerical implementations that beat traditional numerical code, 
either in specialized scientific applications or when compared with popular numerical libraries. The 
following issues are most significant. 

• Awkward e dependence. The sampling complexity, i.e., the number of columns or rows needed 
by the algorithm, scales as or 1/e, which is the usual poor asymptotic complexity for Monte 
Carlo methods. Even though there exist lower bounds for these problems for certain models of data 
access [135j , this is problematic if one wants to choose e to be on the order of machine precision. 

• Numerical conditioning and preconditioning. Thus far, nothing has been said about nu- 
merical conditioning issues, although it is well-known that these issues are crucial when matrix 
algorithms are implemented numerically. 

• Forward error versus backward error. The bounds above, e.g., (fTT|) and (fT2|) . provide so-called 
forward error bounds. The standard stability analysis in NLA is in terms of the backward error, 
where the approximate solution Xopt is shown to be the exact solution of some slightly perturbed 
problem Xopt — argmin^^jKA + SA)x — b\\2, where \ \SA\ \ < e\\A\\ for some small e. 

All three of these issues can be dealt with by coupling the randomized algorithms of the previous sub- 
section with traditional tools from iterative NLA algorithms (as opposed to applying a traditional NLA 
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algorithm as a black box on the random sample or random projection, as the above algorithms do). This 
was first done by and these issues were addressed in much greater detail by [1] and then by [136) 
and [HTI . 

Both of the implementations of [83l [1] take the following form. 

• Premultiply Ahy a structured random projection, e.g., = SHD from Section FS.li which represents 
uniformly sampling a few rows from a randomized Hadamard transform. 

• Perform a QR decomposition on ilA, so that Q.A = QR. 

• Use the R from the QR decomposition as a preconditioner for an iterative Krylov-subspace [55] 
method. 

In general, iterative algorithms compute an e-approximate solution to a LS problem like (jlOp by per- 
forming 0(k(A) log(l/e)) iterations, where k{A) — is the condition number of the input matrix 
(which could be quite large, thereby leading to slow convergence of the iterative algorithm)!^ In this 
case, by choosing the dimension of the random projection appropriately, e.g., as discussed in the previous 
subsections, one can show that k{AR~'-^) is bounded above by a small data-independent constant. That 
is, by using the R matrix from a QR decomposition of VIA, one obtains a good preconditioner for the orig- 
inal problem (|10p . independent of course of any assumptions on the original matrix A. Overall, applying 
the structured random projection in the first step takes o{mii?) time; performing a QR decomposition of 
^A is fast since VIA is much smaller than A; and one needs to perform only 0(log(l/e)) iterations, each 
of which needs 0{mn) time, to compute the approximation. 

The algorithm of [83] used CGLS (Conjugate Gradient Least Squares) as the iterative Krylov-subspace 
method, while the algorithm of [T used LSQR [137'; and both demonstrate that randomized algorithms 
can outperform traditional deterministic NLA algorithms in terms of clock-time (for particular implemen- 
tations or compared with Lapack) for computing high-precision solutions for LS systems with as few as 
thousands of constraints and hundreds of variables. The algorithm of [1] considered five different classes 
of structured random projections (i.e.. Discrete Fourier Transform, Discrete Cosine Transform, Discrete 
Hartely Transform, Walsh-Hadamard Transform, and a Kac random walk), explicitly addressed condi- 
tioning and backward stability issues, and compared their implementation with Lapack on a wide range 
of matrices with uniform, moderately nonuniform, and very nonuniform leverage score importance sam- 
pling distributions. Similar ideas have also been applied to other common NLA tasks; for example, [136] 
shows that these ideas can be used to develop fast randomized algorithms for computing projections onto 
the null space and row space of A, for structured matrices A such that both A and can be rapidly 
applied to arbitrary vectors. 

The implementations of [136[ 1131] are similar, except that the random projection matrix in the first 
step of the above procedure is a traditional Gaussian random projection matrix. While this does not 
lead to a o{mv?) running time, it can be appropriate in certain situations: for example, if both A and its 
adjoint A^ are structured such that they can be applied rapidly to arbitrary vectors [136j : or for solving 
large-scale problems on distributed clusters with high communication cost [ 131j . For example, due to the 
Gaussian random projection, the preconditioning phase of the algorithm of [131] is very well-conditioned, 
which implies that the number of iterations is fully predictable when LSQR or the Chebyshev semi- 
iterative method is applied to the preconditioned system. The latter method is more appropriate for 
parallel computing environments, where communication is a major issue, and thus [131] illustrates the 
empirical behavior of the algorithm on Amazon Elastic Compute Cloud (EC2) clusters. 



These iterative algorithms replace the solution of mot with two problems: first solve Xopt = argmin^||j4n~^?/ — 6||2 
iteratively, where 11 is the preconditioner; and then solve IIx = y. Thus, a matrix fl is a good preconditioner if k{AII^^) 
is small and if Tlx = y can be solved quickly. 
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5 Randomized algorithms for low-rank matrix approximation 

In this section, we will describe several related randomized algorithms for low-rank matrix approximation 
that underlie applications such as those described in Section [2] These algorithms build on the algorithms 
of Section [3.2( and they achieve much- improved worst-case bounds and are more useful in both numerical 
analysis and data analysis applications (when compared with the algorithms of Section . First, the 
algorithm of Section fS.ll is a random sampling algorithm that improves the additive-error bounds to much 
finer relative-error bounds — the analysis of this algorithm boils down to an immediate application of the 
least-squares approximation algorithm of Section |4l Next, the algorithm of Section [5.2l is a random sam- 
pling algorithm that returns exactly k, rather than 0(fclogfc/e^), columns, for an input rank parameter 
k — the proof of this result involves a structural property that decouples the randomization from the linear 
algebra in a somewhat more refined way than we saw in Section Finally, the algorithms of Section 15.31 
are random projection algorithms that take advantage of this more refined structural property to couple 
these randomized algorithms with traditional methods from NLA and scientific computingo 



5.1 A basic random sampling algorithm 

Additive-error bounds (of the form proved for the low-rank algorithms of Section 13. 2p are rather coarse, 
and the gold standard in TCS is to establish much finer relative-error bounds of the form provided 
in (jl6p below. To motivate the importance sampling probabilities used to achieve such relative-error 
guarantees, recall that if one is considering a matrix with fc — 1 large singular values and one much 
smaller singular value, then the directional information of the fc*'' singular direction will be hidden from 
the Euclidean norms of the input matrix. The reason is that, since Ak = Uk^kyk , the Euclidean norms 
of the columns of A are convolutions of "subspace information" (encoded in Uk and V^) and "size- 
oi-A information" (encoded in Sfe). This suggests deconvoluting subspace information and size-of-A 
information by choosing importance sampling probabilities that depend on the Euclidean norms of the 
columns of . This importance sampling distribution defines a non- uniformity structure over R" that 
indicates where in the n-dimensional space the information in A is being sent, independent of what that 
(singular value) information is. More formally, these quantities are proportional to the diagonal elements 
of the projection matrix onto the span of V^"'"!^^ and thus they are examples of generalized statistical 
leverage scores. 

This idea was suggested in [1151 [2T| . and it forms the basis for the algorithm from the TCS literature 
that achieves the strongest Frobenius norm boundsl^^ Given an m x n matrix A and a rank parameter fc: 

• Compute the importance sampling probabilities where pi — ^||V'^"^'''*||2, where VjJ is any 
k X n orthogonal matrix spanning the top-fc right singular subspace of A. 

• Randomly select and rescale c = 0{k log fc/e^) columns of A according to these probabilities to form 



■^^Most of the results in this section will be formulated in terms of the amount of spectral or Frobenius norm that is 
captured by the (full or rank-fc approximation to the) random sample or random projection. Given a basis for this sample 
or projection, it is straightforward to compute other common decompositions such as the pivoted QR factorization, the 
eigenvalue decomposition, the partial SVD, etc. using traditional NLA methods; see |14| for a good discussion of this. 

'^''Here, we are sampling columns and not rows, as in the algorithms of Section |4l and thus we are dealing with the right, 
rather than the left, singular subspace; but clearly the ideas are analogous. Thus, in particular, note that the "span of 
yj" refers to the span of the rows of V^, whereas the "span of Uk," as used in previous sections, refers to the span of the 
columns of 11^. 

Subsequent work in TCS that has not (yet?) found application in NLA and data analysis also achieves similar relative- 
error bounds but with different methods. For example, the algorithm of 138 runs in roughly 0(mnk^ log k) time and uses 
geometric ideas involving sampling and merging approximately optimal fc-flats. Similarly, the algorithm of [1391 randomly 
samples in a more complicated manner and runs in 0{Mk^ log fc), where M is the number of nonzero elements in the matrix; 
alternatively, it runs in O(fclogfc) passes over the data from external storage. 
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the matrix C . 

A more detailed description of this basic random samphng algorithm may be found in |115[ [2T| , where it 
is proven that 

U - Pc,A\\p < (1 + e)\\A - Pu,A\\f (16) 

holds. (As above, Ck is the best rank-/c approximation to the matrix C, and Pc^. is the projection matrix 
onto this fc-dimensional space.) As with the relative-error random sampling algorithm of Section [4. 2 1 the 
dependence of the sampling complexity and running time on the failure probability 6 is 0(log(l/(5)); thus, 
the failure probability for this randomized low-rank approximation, as well as the subsequent algorithms 
of this section, can be made to be negligibly-small, both in theory and in practice. The analysis of this 
algorithm boils down to choosing a set of columns that are relative-error good at capturing the Frobenius 
norm of A, when compared to the basis provided by the top-fc singular vectors. That is, it boils down 
to the randomized algorithm for the least-squares approximation problem from Section |4l see [115| 121) 
for details. 

This algorithm and related algorithms that randomly sample columns and/or rows provide what is 
known as CX or CUR matrix decompositions [1141 11151 121) rn In addition, this relative-error column- 
based CUR decomposition, as well as heuristic variants of it, has been applied in a range of data analysis 
applications, ranging from term-document data to DNA SNP data [21] [25l [10] . The computational bot- 
tleneck for this relative-error approximation algorithm is computing the importance sampling distribution 
{Pi}?=ii for which it suffices to compute any kxn matrix that spans the top-A; right singular subspace 
of A. That is, it suffices (but is not necessary) to compute any orthonormal basis spanning , which 
typically requires 0{mnk) running time, and it is not necessary to compute the full or partial SVD. 
Alternatively, the leverage scores can all be approximated to within 1 ± e in roughly 0{mnk) time using 
the algorithm of 133J from Section [4.4.21 in which case these approximations can be used as importance 
sampling probabilities in the above random sampling algorithm. 

5.2 A more refined random sampling algorithm 

The algorithm of the previous subsection randomly samples 0{k\ogk/e^) columns, and then in order 
to compare with the bound provided by the SVD, it "filters" those columns through an exactly rank-/c 
space. In this subsection, we will describe a randomized algorithm for the Column Subset Selection 
Problem (CSSP): the problem of selecting exactly k columns from an input matrix. Clearly, bounds 
for algorithms for the CSSP will be worse than the very strong relative-error bounds, provided by ([T6|) . 
that hold when 0(A:logfc/e^) columns are selected. Importantly, though, this CSSP algorithm extends 
the ideas of the previous relative-error TCS algorithm to obtain bounds of the form used historically in 
NLA. Moreover, the main structural result underlying the analysis of this algorithm permits much finer 
control on the application of randomization, e.g., for high-performance numerical implementation of both 
random sampling and random projection algorithms. We will start with a review of prior approaches to 
the CSSP; then describe the algorithm and the quality-of-approximation bounds; and finally highlight 
the main structural result underlying its analysis. 

5.2.1 A formalization of and prior approaches to this problem 

Within NLA, a great deal of work has focused on this CSSP [Mil [Hg [1461 HHl [UHl [HSl HlOl [IS [TKTl 
I152|, 1153] . Several general observations about the NLA approach include: 



•^^Within the NLA community, Stewart developed the quasi-Gram-Schmidt method and apphed it to a matrix and its 
transpose to obtain such a CUR matrix decomposition |140lll4l| : and Goreinov, Tyrtyshnikov, and Zamarashkin developed 
a CUR matrix decomposition (a so-called pseudoskeleton approximation) and related the choice of columns and rows to a 
"maximum uncorrelatedness" concept |142II143] . Note that the Nystrom method is a type of CUR decomposition and that 
the pseudoskeleton approximation is also a generalization of the Nystrom method. 
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• The focus in NLA is on deterministic algorithms. Moreover, these algorithms are greedy, in that 
at each iterative step, the algorithm makes a decision about which columns to keep according to 
a pivot-rule that depends on the columns it currently has, the spectrum of those columns, etc. 
Differences between different algorithms often boil down to how to deal with such pivot rules 
decisions, and the hope is that more sophisticated pivot-rule decisions lead to better algorithms in 
theory or in practice. 

• There are deep connections with QR factorizations and in particular with the so-called Rank Reveal- 
ing QR factorizations. Moreover, there is an emphasis on optimal conditioning questions, backward 
error analysis issues, and whether the running time is a large or small constant multiplied by 

or n^. 

• Good spectral norm bounds are obtained. A typical spectral norm bound is: 

\\A~PcA\\2<0(^^k{n-k)) \\A-Pu,A\\2, (17) 

and these results are algorithmic, in that the running time is a low-degree polynomial in m and 
n [13j . (The ^/n multiplicative factor might seem large, but recall that the spectral norm is much 
less "forgiving" than the Frobenius norm and that it is not even known whether there exists columns 
that do better.) On the other hand, the strongest result for the Frobenius norm in this literature is 

\\A - PcA\\f < V{k + l){n-k)\\A ~ Pu,A\\2, (18) 

but it is only an existential result, i.e., the only known algorithm essentially involves exhaustive 
enumeration [149j . (In these two expressions, Uk is the m x k matrix consisting of the top k left 
singular vectors of A, and Pjj^. is a projection matrix onto the span of Uk-) 

Within TCS, a great deal of work has focused on the related problem of choosing good columns from 
a matrix [921 [60l [94l 11141 l65l 1115) . Several general observations about the TCS approach include: 

• The focus in TCS is on randomized algorithms. In particular, with these algorithms, there exists 
some nonzero probability, which can typically be made extremely small, say 6 — 10~^°, that the 
algorithm will return columns that fail to satisfy the desired quality-of-approximation bound. 

• The algorithms select more than k columns, and the best rank-fc projection onto those columns is 
considered. The number of columns is typically a low-degree polynomial in k, most often 0{k log k), 
where the constants hidden in the big-0 notation are quite reasonable. 

• Very good Frobenius norm bounds are obtained. For example, the algorithm (described above) that 
provides the strongest Frobenius norm bound achieves (jl6p . while running in time of the order of 
computing an exact or approximate basis for the top-fc right singular subspace |115j . The TCS 
literature also demonstrates that there exists a set of k columns that achieves a constant-factor 
approximation: 

\\A-PcA\\F<Vk\\A^Pu,A\\F, (19) 
but note that this is an existential result |154j . 

Note that, prior to the algorithm of the next subsection, it was not immediately clear how to combine 
these two very different approaches. For example, if one looks at the details of the pivot rules in the 
deterministic NLA methods, it isn't clear that keeping more than exactly k columns will help at all 
in terms of reconstruction error. Similarly, since there is a version of the so-called Coupon Collector 
Problem at the heart of the usual TCS analysis, keeping fewer than 57(fclogfc) will fail with respect to 
this worst-case analysis. Moreover, the obvious hybrid algorithm of first randomly sampling O(fclogfc) 
columns and then using a deterministic QR procedure to select exactly k of those columns does not seem 
to perform so well (either in theory or in practice). 
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5.2.2 A two-stage hybrid algorithm for this problem 

Consider the following more sophisticated version of a two-stage hybrid algorithm. Given an arbitrary 
m X n matrix A and rank parameter k: 

• (Randomized phase) Compute the importance sampling probabilities {pi}^^i, where = -^H V^"'"*'*'' Hi; 
where VlJ" is any kxn orthogonal matrix spanning the top-fc right singular subspace of A. Randomly 
select and rescale c = 0{k log k) columns of Vj[ according to these probabilities. 

• (Deterministic phase) Let V"^ be the k x O(fclogfc) non-orthogonal matrix consisting of the down- 
sampled and rescaled columns of V^f . Run a deterministic QR algorithm on V'^ to select exactly k 
columns of . Return the corresponding columns of A. 

A more detailed description of this algorithm r nay be found in [24l I116| . where it is shown that with 
extremely high probability the following spectrao and Frobenius norm bounds hold: 

\\A-PcA\\2 < Oik'/Hog'/\ky^')\\A-Pu,A\\2 (20) 
\\A~PcA\\f < 0{klog'/H)\\A~Pu,A\\F. (21) 

Note that both the original choice of columns in the first phase, as well as the application of the QR 
algorithm in the second phased involve the sampled version of the matrix Vjf , i.e., the matrix defining 
the relevant non- uniformity structure over the columns of A in the relative-error algorithm of Section [5. II 
In particular, it is critical to the success of this algorithm that the QR procedure in the second phase be 
applied to the randomly-sampled version of Vj[ , rather than of A itself. This algorithm may be viewed 
as post-processing the relative-error random sampling algorithm of the previous subsection to remove 
redundant columns; and it has been applied successfully to a range of data analysis problems. See, 
e.g., [inillllllS] and |155l 11561 1157| . as well as |158j for a discussion of numerical issues associated with 
this algorithm. 

With respect to running time, the computational bottleneck for this algorithm is computing {pi}f^i, 
for which it suffices to compute any kxn matrix that spans the top-k right singular subspace of 
A. (In particular, a full or partial SVD computation is not necessary.) Thus, this running time is 
of the same order as the running time of the QR algorithm used in the second phase when applied 
to the original matrix A, typically roughly 0{mnk) time. (Not surprisingly, one could also perform a 
random projection, such as those described in Section 15.31 below, to approximate this basis, and then 
use that approximate basis to compute approximate importance sampling probabilities, as described in 
Section [4.4. 21 above. In that case, similar bounds would hold, but the running time would be improved to 
0(mn log fc) time.) Moreover, this algorithm scales up to matrices with thousands of rows and millions 
of columns, whereas existing off-the-shelf implementations of traditional QR algorithms may fail to run 
at all. With respect to the worst-case quality of approximation bounds, this algorithm selects columns 
that are comparable to the state-of-the-art algorithms for constant k {i.e., ((20|) is only 0(/c^/* log^^^ k) 
worse than (fT7|) ) for the spectral norm; and (j21|l is only a factor of at most 0((A; log /c)^/^) worse than 
(jl9p . the best previously- known existential result for the Frobenius norm. 

5.2.3 A basic structural result 

As with the relative-error LS algorithm of Section |4l in order to see "why" this algorithm for the CSSP 
works, it is helpful to identify a structural condition that decouples the randomization from the linear 



•^•^Note that to establish the spectral norm bound, |24l 11161 used slightly more complicated (but still depending only on 
information in V^) importance sampling probabilities, but this may be an artifact of the analysis. 

^*Note that QR (as opposed to the SVD) is not performed in the second phase to speed up the computation of a relatively 
cheap part of the algorithm, but instead it is performed since the goal of the algorithm is to return actual columns of the 
input matrix. 
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algebra. This structural condition was first identified in [531 1116) . and it was subsequently improved 
by |14| . To identify it, consider preconditioning or postmultiplying the input matrix A by some arbitrary 
matrix Z. Thus, for the above randomized algorithm, the matrix Z is a carefully-constructed random 
sampling matrix, but it could be a random projection, or more generally any other arbitrary matrix Z . 
Recall that if fc < r = rank(A), then the SVD of A may be written as 



where Uk is the m x k matrix consisting of the top k singular vectors, Uk,± is the m x (r — k) matrix 
consisting of the bottom r — k singular vectors, etc. Then, the following structural condition holds. 

• Structural condition underlying the randomized low-rank algorithm. If V^f Z has full 
rank, then for i/ G {2, F}, i.e., for both the Frobenius and spectral norms. 



holds, where Paz is a projection onto the span of AZ , and where the dagger symbol represents the 
Moore-Penrose pseudoinverse. 

This structural condition characterizes the manner in which the behavior of the low-rank algorithm 
depends on the interaction between the right singular vectors of the input matrix and the matrix Z. 
(In particular, it depends on the interaction between the subspace associated with the top part of the 

spectrum and the subspace associated with the bottom part of the spectrum via the (vi^^^Z^ {V^Zy 

term.) Note that the assumption that VjF Z does not lose rank is a generalization of Condition of 
Section 31 Also, note that the form of this structural condition is the same for both the spectral and 
Frobenius norms. 

As with the LS problem, given this structural insight, what one does with it depends on the application: 
one can compute the basis Vjf exactly if that is not computationally prohibitive and if one is interested 
in extracting exactly k columns; or one can perform a random projection and ensure that with high 
probability the structural condition is satisfied. Moreover, by decoupling the randomization from the 
linear algebra, it is easier to parameterize the problem in terms more familiar to NLA and scientific 
computing: for example, one can consider sampling £ > k columns and projecting onto a rank-A;', where 
fc' > k, approximation to those columns; or one can couple these ideas with traditional methods such as 
the power iteration method. Several of these extensions will be the topic of the next subsection. 

5.3 Several related random projection algorithms 

In this subsection, we will describe three random projection algorithms that draw on the ideas of the 
previous subsections in progressively finer ways. 

5.3.1 A basic random projection algorithm 

To start, consider the following basic random projection algorithm. Given an m x n matrix A and a rank 
parameter fc: 

• Construct an n x £, with £ = 0{k/e), structured random projection matrix f2, e.g., fl = DHS 
from Section 13. li which represents uniformly sampling a few rows from a randomized Hadamard 
transform. 

• Return B = AO. 




(22) 
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This algorithm, which amounts to choosing uniformly a small number £ of columns in a randomly rotated 
basis, was introduced in [80], where it is proven that 

\\A-PB,A\\F<{l + e)\\A-Pu,A\\F (23) 

holds with high probability. (Recall that Bk is the best rank-fc approximation to the matrix B, and PB^. 
is the projection matrix onto this /c-dimensional space.) This algorithm runs in 0{Mk/e+ (m + n)k'^ / e^) 
time, where M is the number of nonzero elements in A, and it requires 2 passes over the data from 
external storage. 

Although this algorithm is very similar to the additive-error random projection algorithm of |91| that 
was described in Section 13. 2[ this algorithm achieves much stronger relative-error bounds by performing 
a much more refined analysis. Basically, [80] (and also the improvement [159j ) modifies the analysis of 
the relative-error random sampling of [115| 121) that was described in Section 15.11 which in turn relies 
on the relative-error random sampling algorithm for LS approximation [1271 1115j that was described in 
Sectional In the same way that we saw in Section [4.41 that fast structured random projections could be 
used to uniformize coordinate-based non-uniformity structure for the LS problem, after which fast uniform 
sampling was appropriate, here uniform sampling in the randomly-rotated basis achieves relative-error 
bounds. In showing this, [80] also states a "subspace" analogue to the JL lemma, in which the geometry 
of an entire subspace of vectors (rather than just N pairs of vectors) is preserved. Thus, one can view 
the analysis of [50] as applying JL ideas, not to the rows of A itself, as was done by [HI], but instead to 
vectors defining the subspace structure of A. Thus, with this random projection algorithm, the subspace 
information and sizc-of-A information are deconvoluted within the analysis, whereas with the random 
sampling algorithm of Section 15.11 this took place within the algorithm by modifying the importance 
sampling probabilities. 

5.3.2 An improved random projection algorithm 

As with the randomized algorithms for the LS problem, several rubber-hits-the-road issues need to be 
dealt with in order for randomized algorithms for the low-rank matrix approximation problem to yield to 
high-precision numerical implementations that beat traditional deterministic numerical code. In addition 
to the issues described in Section [45] the main issue here is the following. 

• Minimizing the oversampling factor. In practice, choosing even 0{k\ogk) columns, in either 
the original or a randomly-rotated basis, even when the big-0 notation hides only modest constants, 
can make it difhcult for these randomized matrix algorithms to beat previously-existing high-quality 
numerical implementations. Ideally, one could parameterize the problem so as to choose some 
number i = k + p columns, where p is a modest additive oversampling factor, e.g., 10 or 20 or k, 
and where there is no big-0 constant. 

When attempting to be this aggressive at minimizing the size of the sample, the choice of the oversampling 
factor p is more sensitive to the input than in the algorithms we have reviewed so far. That is, whereas 
the previous bounds held for any input, here the proper choice for the oversampling factor p can be 
quite sensitive to the input matrix. For example, when parameterized this way, p could depend on the 
size of the matrix dimensions, the decay properties of the spectrum, and the particular choice made 
for the random projection matrix [1601 11611 [5^ 11621 HM [TB5] . Moreover, for worst-case input matrices, 
such a procedure may fail. For example, one can very-easily construct matrices such that if one randomly 
samples o(k log k) columns, in either the original canonical basis or in the randomly-rotated basis provided 
by the structured Hadamard transform, then the algorithm will fail. Basically, one can easily encode the 
so-called Coupon Collector Problem [128j in the columns, and it is known that O(fclogfc) samples are 
necessary for this problem. 

That being said, running the risk of such a failure might be acceptable if one can efficiently couple 
to a diagnostic to check for such a failure, and if one can then correct for it by choosing more samples 
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if necessary. The best numerical implementations of randomized matrix algorithms for low-rank matrix 
approximation do just this, and the strongest results in terms of minimizing p take advantage of Con- 
dition (1^^ in a somewhat different way than was originally used in the analysis of the CSSP [13]. For 
example, rather than choosing 0{k log k) dimensions and then filtering them through exactly k dimen- 
sions, as the relative-error random sampling and relative-error random projection algorithms do, one can 
choose some number i of dimensions and project onto a fc'-dimensional subspace, where k < k' < i, 
while exploiting Condition (1221) to bound the error, as appropriate for the computational environment at 
hand [M]. 

Next, consider a second random projection algorithm that will address this issue. Given an m x n 
matrix A, a rank parameter fc, and an oversampling factor p: 

• Set £ = k+p. 

• Construct an n x £ random projection matrix $7, either with i.i.d. Gaussian entries or in the form 
of a structured random projection such as = DHS which represents uniformly sampling a few 
rows from a randomized Hadamard transform. 

• Return B ^ AVt 

Although this is quite similar to the algorithms of |91[ 180] , algorithms parameterized in this form were 
introduced in |160| 11611 [82] , where a suite of bounds of the form 



are shown to hold with high probability. Here, Z is a rank-Zc-or-greater matrix easily-constructed from B. 
This result can be used to obtain a so-called interpolative decomposition (a variant of the basic CSSP with 
explicit numerical conditioning properties) , and jl60[ 11611 15^ also provide an a posteriori error estimate 
(that is useful for situations in which one wants to choose the rank parameter k to be the numerical rank, 
as opposed to the a priori specification of k as part of the input, which is more common in the TCS-style 
algorithms that preceded this algorithm). 

5.3.3 A third random projection algorithm 

Finally, consider a third random projection algorithm that will address the issue that the decay properties 
of the spectrum can be important when it is of interest to minimize the oversampling very aggressivelvl^ 
Given an m x n matrix A, a rank parameter k, an oversampling factor p, and an iteration parameter q: 

• Set £ = k+p. 

• Construct an n x ^ random projection matrix Q, either with i.i.d. Gaussian entries or in the form 
of a structured random projection such as f2 = DHS which represents uniformly sampling a few 
rows from a randomized Hadamard transform. 

• Return B = {AA^)''An 

This algorithm (as well as a numerically-stable variant of it) was introduced in [162| . where it is shown 
that bounds of the form 



^^Of course, this should not be completely unexpected, given that Condition I I22II shows that the behavior of algorithms 
depends on the interaction between different subspaces associated with the input matrix A. When stronger assumptions 
are made about the data, stronger bounds can often be obtained. 



\\A- Z\\2 < 10y^£mm{m,n}\\A- A, 



-k\\2 
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hold with high probabiUty. (This bound should be compared with the bound for the previous algorithm, 
and thus Z is a rank-fc-or-greater matrix easily-constructed from B.) Basically, this random projection 
algorithm modifies the previous algorithm by coupling a form of the power iteration method within the 
random projection step. This has the effect of speeding up the decay of the spectrum while leaving the 
singular vectors unchanged, and it is observed in |162| [T4] that q = 2 or g = 4 is often sufficient for 
certain data matrices of interest. This algorithm was analyzed in greater detail for the case of Gaussian 
random matrices in |14) , and an out-of-core implementation (meaning, appropriate for data sets that are 
too large to be stored in RAM) of it was presented in [163] . 

The running time of these last two random projection algorithms depends on the details of the 
computational environment, e.g., whether the matrix is large and dense but fits into RAM or is large 
and sparse or is too large to fit into RAM; how precisely the random projection matrix is constructed; 
whether the random projection is being applied to an arbitrary matrix A or to structured input matrices, 
etc. |14| . For example, if random projection matrix is constructed from i.i.d. Gaussian entries then 
in general the algorithm requires 0{mnk) time to implement the random projection, i.e., to perform 
the matrix-matrix multiplication AQ, which is no faster than traditional deterministic methods. On 
the other hand, if the projection is to be applied to matrices A such that A and/or A'^ can be applied 
rapidly to arbitrary vectors {e.g., very sparse matrices, or structured matrices such as those arising from 
Toeplitz operators, or matrices that arise from discretized integral operators that can be applied via the 
fast multipole method), then Gaussian random projections may be preferable. Similarly, in general, if fl 
is structured, e.g., is of the form fl = DHS, then it can be implemented in O(TOnlogfc) time, and this 
can lead to dramatic clock-time speed-up over classical techniques even for problems of moderate sizes. 
On the other hand, for out-of-core implementations these additional speed-ups have a negligible effect, 
e.g., since matrix- matrix multiplications can be faster than a QR factorization, and so using Gaussian 
projections can be preferable. Working through these issues in theory and practice is still very much an 
active research area. 



6 Empirical observations 

In this section, we will make some empirical observations, with an emphasis on the role of statistical 
leverage in these algorithms and in MMDS applications more generally. 



6.1 Traditional perspectives on statistical leverage 

As mentioned previously, the use of statistical leverage scores has a long history in statistical and diag- 
nostic regression analysis [1191 11181 1120[ I12H I122j . To gain insight into these statistical leverage scores. 



consider the so-called "wood beam data" example [1641 1119[ , which is visually presented in Figure 2(a) 
along with the best-fit line to that data. In Figure [2(b)[ the leverage scores for these ten data points are 
shown. Intuitively, data points that "stick out" have particularly high leverage — e.g., the data point that 
has the most influence or leverage on the best-fit line to the wood beam data is the point marked "4" , and 
this is reflected in the relative magnitude of the corresponding statistical leverage score. (Note that the 
point "1" exhibits some similar behavior; and that although "3" and "9" don't "stick out" in the same 
sense, they are at the "ends" of the data and possess a relatively-high leverage for that reason.) Indeed, 
since Trace(i?) = n, where H is the hat matrix defined in Section 14. li a rule of thumb that has been 
suggested in diagnostic regression analysis to identify errors and outliers in a data set is to investigate 
the data point if Hu > 2n/m \121\ 1122] . i.e., if Ha is larger that 2 or 3 times the "average" size. On 
the other hand, of course, if it happens to turn out that such a point is a legitimate data point, then one 
might expect that such an outlying data point will be a particularly important or informative data point. 

That leverage scores "should be" fairly uniform — indeed, typical conditions even in recent work on the 
so-called coherence of matrices [165i 11661 HI 1125] make just such an assumption — is supported by the idea 
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Figure 2: Statistical leverage scores historically in diagnostic regression analysis. (2(a) ) The Wood Beam 
Data described in |119j is an example illustrating the use of statistical leverage scores in the context 
of least-squares approximation. Shown are the original data and the best least-squares fit. (2(b) I The 
leverage scores for each of the ten data points in the Wood Beam Data. Note that the data point marked 
"4" has the largest leverage score, as might be expected from visual inspection. 



that if they are not then a small number of data points might be particularly important, in which case a 
different or more refined statistical model might be appropriate. Furthermore, they are fairly uniform in 
various limiting cases where measure concentration occurs, e.g., for not-extremely-sparse random graphs, 
and for matrices such as Laplacians associated with well-shaped low-dimensional manifolds, basically since 
eigenfunctions tend to be delocalized in those situations. Of course, their actual behavior in realistic data 
applications is an empirical question. 



6.2 More recent perspectives on statistical leverage 

To gain intuition for the behavior of the statistical leverage scores in a typical application, consider 



Figure 3(a) which illustrates the so-called Zachary karate club network 167 , a small but popular network 
in the community detection literature. Given such a network G = {V,E), with n nodes, m edges, and 
corresponding edge weights We > define the n x n Laplacian matrix as i = B^WB, where B is 
the m X n edge-incidence matrix and W is the m x m diagonal weight matrix. The effective resistance 
between two vertices is given by the diagonal entries of the matrix R = BL^ B^ (where denotes the 
Moore-Penrose generalized inverse) and is related to notions of "network betweenness" |168| . For many 
large graphs, this and related betweenness measures tend to be strongly correlated with node degree and 
tend to be large for edges that form articulation points between clusters and communities, i.e., for edges 
that "stick out" a lot. It can be shown that the effective resistances of the edges of G are proportional 
to the statistical leverage scores of the m rows of the mx n matrix W^/'^B — consider the mx m matrix 

P = W^/'^RW^''^ = $($^$)^$'^, 

where $ = W^^^B, and note that if C/$ denotes any orthogonal matrix spanning the column space of 
then 

Pa — iU<s>U^)ii = ||(f^*)(i)|l2- 



Figure 3(a) presents a color-coded illustration of these scores for Zachary karate club network. Note 



that the higher-leverage red edges tend to be those associated with higher-degree nodes and those at the 
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"Cumulative leverage" for terms In "Enron Matrix" 




(a) Zachary Karate Club Data (b) Cumulative Leverage 
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(c) Leverage Score and Information Gain for DNA Microarray Data 



Figure 3: Statistical leverage scores in more modern applications. (3(a)) The so-called Zachary karate 
club network [167^, with edges color-coded such that leverage scores for a given edge increase from 
yellow to red. (3(b) I Cumulative leverage (with k = 10) for a 65,031 x 92,133 term-document matrix 
constructed Enron electronic mail collection, illustrating that there are a large number of data points 
with very high leverage score. (3(c) I The normalized statistical leverage scores and information gain 
score — information gain is a mutual information-based metric popular in the application area [lOi 121] — 
for each of the n = 5520 genes, a situation in which the data cluster well in the low-dimensional space 
defined by the maximum variance axes of the data jH]. Red stars indicate the 12 genes with the highest 
leverage scores, and the red dashed line indicates the average or uniform leverage scores. Note the strong 
correlation between the unsupervised leverage score metric and the supervised information gain metric. 



articulation point between the two clusters. 

Next, to gain intuition for the (non-) uniformity properties of statistical leverage scores in a typical 
application, consider a term-document matrix derived from the publicly-released Enron electronic mail 
collection [169] . which is typical of the type of data set to which SVD-based latent semantic analysis 
(LSA) methods |170j have been applied. This is a 65,031 x 92, 133 matrix, as described in [169|, and let 
us choose the rank parameter as fc = 10. Figure [3 (b) | plots the cumulative leverage, i.e., the running sum 
of top t statistical leverage scores, as a function of increasing t. Since ^ = g^j^ ~ 1.0854 x 10""*, we 
see that the highest leverage term has a leverage score nearly two orders of magnitude larger than this 
"average" size scale, that the second highest-leverage score is only marginally less than the first, that the 
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third highest score is marginahy less than the second, etc. Thus, by the traditional metrics of diagnostic 
data analysis |121i 1122] , which suggests flagging a data point if 

there are a huge number of data points that are extremely outlying, i.e., that are extreme outliers by the 
metrics of traditional regression diagnostics. In retrospect, of course, this might not be surprising since 
the Enron email corpus is extremely sparse, with nowhere on the order of n{n) nonzeros per row. Thus, 
even though LSA methods have been successfully applied, plausible generative models associated with 
these data are clearly not Gaussian, and the sparsity structure is such that there is no reason to expect 
that nice phenomena such as measure concentration occur. 

Finally, note that DNA microarray and DNA SNP data often exhibit a similar degree of non- 



uniformity, although for somewhat different reasons. To illustrate. Figure 3(c) presents two plots for 
a data matrix, as was described in consisting of m = 31 patients with 3 different cancer types with 
respect to n = 5520 genes. First, this figure plots the information gain — information gain is a mutual 
information-based metric popular in that application area [TOl [H] — for each of the n = 5520 genes; and 
second, it plots the normalized statistical leverage scores for each of these genes. In each case, red dots 
indicate the genes with the highest values. A similar plot illustrating the remarkable non-uniformity in 
statistical leverage scores for DNA SNP data was presented in [TIT . Empirical evidence suggests that two 
phenomena may be responsible for this non-uniformity. First, as with the term-document data, there is 
no domain-specific reason to believe that nice properties like measure concentration occur — on the con- 
trary, there are reasons to expect that they do not. Recall that each DNA SNP corresponds to a single 
mutational event in human history. Thus, it will "stick out," as its description along its one axis in the 
vector space will likely not be well-expressed in terms of the other axes, i.e., in terms of the other SNPs, 
and by the time it "works its way back" due to population admixing, etc., other SNPs will have occurred 
elsewhere. Second, the correlation between statistical leverage and supervised mutual information-based 
metrics is particularly prominent in examples where the data cluster well in the low- dimensional space 
defined by the maximum variance axes. Considering such data sets is, of course, a strong selection bias, 
but it is common in applications. It would be of interest to develop a model that quantifies the obser- 
vation that, conditioned on clustering well in the low-dimensional space, an unsupervised measure like 
leverage scores should be expected to correlate well with a supervised measure like informativeness [lO] 
or information gain |21j . 



6.3 Statistical leverage and selecting columns from a matrix 

With respect to some of the more technical and implementational issues underlying the CSSP algorithm 
of Section [SI recall that an important aspect of QR algorithms is how they make so-called pivot rule 
decisions about which columns to keep |53) and that such decisions can be tricky when the columns 
are not orthogonal or spread out in similarly nice ways. Several empirical observations |155[ 1116) are 
particularly relevant for large-scale data applications. 

• We looked at several versions of the QR algorithm, and we compared each version of QR to the 
CSSP using that version of QR in the second phase. One observation we made was that different 
QR algorithms behave differently — e.g., some versions such as the Low-RRQR algorithm of (171) 
tend to perform much better than other versions such as the qrxp algorithm of |151i 1172) . Although 
not surprising to NLA practitioners, this observation indicates that some care should be paid to 
using "off the shelf" implementations in large-scale applications. A second less-obvious observation 
is that preprocessing with the randomized first phase tends to improve more poorly-performing 
variants of QR more than better variants. Part of this is simply that the more poorly-performing 
variants have more room to improve, but part of this is also that more sophisticated versions of QR 
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tend to make more sophisticated pivot rule decisions, which are relatively less important after the 
randomized bias toward directions that are "spread out." 

• We also looked at selecting columns by applying QR on and then keeping the corresponding 
columns of A, i.e., just running the classical deterministic QR algorithm with no randomized first 
phase on the matrix VjF . Interestingly, with this "preprocessing" we tended to get better columns 
than if we ran QR on the original matrix A. Again, the interpretation seems to be that, since the 
norms of the columns of V^F define the relevant non-uniformity structure with which to sample with 
respect to, working directly with those columns tends make things "spread out," thereby avoiding 
(even in traditional deterministic settings) situations where pivot rules have problems. 

• Of course, we also observed that randomization further improves the results, assuming that care is 
taken in choosing the rank parameter k and the sampling parameter c. In practice, the choice of k 
should be viewed as a "model selection" question. Then, by choosing c = k, 1.5k, 2k, . . ., we often 
observed a "sweet spot," in bias-variance sense, as a function of increasing c. That is, for a fixed 
k, the behavior of the deterministic QR algorithms improves by choosing somewhat more than k 
columns, but that improvement is degraded by choosing too many columns in the randomized phase. 

These and related observations |1551 1116) shed light on the inner workings of the CSSP algorithm, 
the effect of providing a randomized bias toward high-leverage data points at the two stages of the 
algorithm, and potential directions for the usefulness of this type of randomized algorithm in very large- 
scale data applications. 

6.4 Statistical leverage in large-scale data analysis 

Returning to the genetics applications where the algorithms described in this review have already been 
applied, we will consider one example each of the two common reasons (faster algorithms and more- 
interpretable algorithms) described in Section [^751 for using randomization in the design of matrix algo- 
rithms for large-scale data problems. To start with the former motivation, [T^ applies the algorithm 
of [162] that was described in Section 15.3.31 to problems of subspace estimation and prediction in case- 
control studies in genetics. In particular, [15] extends Principal Component Regression, Sliced Inverse 
Regression, and Localized Sliced Inverse Regression, three statistical techniques for matrix-based data 
that use as a critical step a matrix eigendecomposition, to much larger-scale data by using randomization 
to compute an approximately-optimal basis. Three goals were of interest: evaluate the ability of random- 
ized matrix algorithms to provide a good approximation to the dimension-reduced subspace; evaluate 
their relative predictive performance, when compared to the exact methods; and evaluate their ability to 
provide useful graphical summaries for the domain experts. 

As an example of their results. Figure |4] illustrates pictorially the ability of the randomized algorithm 
to recover a good approximation to the top principal components of the DNA SNP data of ,175] . For 
appropriate parameter choices, the empirical correlations of the exact sample principal components with 
the approximate estimates are (starting from second principal component): 0.999, —0.996, 0.930, 0.523, 
0.420, 0.255, etc. Thus, the first few components are extremely-well reproduced; and then it appears as if 
there is some "component splitting," as would be expected from standard matrix perturbation theory, as 
the seventh principal component from the randomized algorithm correlates relatively-well with three of 
the exact components. The performance of the randomized algorithm in a Crohn's disease application on 
a subset of the DNA SNP data from the Wellcome Trust Case Control Consortium |174j illustrates the 
run-time advantages of exploiting randomization in a real case-control application. For example, when 
applied to a 4,686 x 6,041 matrix of SNPs from one chromosome, a single iteration of the randomized 
algorithm took 6 seconds, versus 37 minutes for the exact deterministic computation (a call to the 
DGESVD routine in the Lapack package), and it achieved a distance of 0.01 to the exact subspace. By 
iterating just a few additional times, this distance could be decreased to less than 10"^ with a nominal 
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Figure 4: Pictorial illustration of the method of [15] recovering a good approximation to the top principal 
components of the DNA SNP data of [173) . Left panel is with the randomized algorithm; and right panel 
is with an exact computation. See the text and 1151 for details. 



increase in running time relative to the exact computation. Similar results were achieved with matrices 
of larger sizes, up to a 4, 686 x 29, 406 matrix consisting of SNP data from four chromosomes. 

In a somewhat different genetics application, |10| was interested in obtaining a small number of actual 
DNA SNPs that could be typed and used for ancestry inference and the study of population structure 
within and across continents around the world. As an example of their results, Figure [5] illustrates 
pictorially the clustering of individuals from nine indigenous populations typed for 9, 160 SNPs. fc-means 
clustering was run on the detected significant eigenvectors, and it managed successfully to assign each 
individual to their country of origin. In order to determine the possibility of identifying a small set of 
actual SNPs to reproduce this clustering structure, [10] used the statistical leverage scores as a ranking 
function and selected sets of 10 to 400 actual SNPs and repeated the analysis. Figure [5] also illustrates 
the correlation coefficient between the true and predicted membership of the individuals in the nine 
populations, when the SNPs are chosen in this manner, as well as when the SNPs are chosen with two other 
procedures (uniformly at random and according to a mutual information-based measure popular in the 
field). Surprisingly, by using only 50 such "PCA-correlated" SNPs, all individuals were correctly assigned 
to one of the nine studied populations, as illustrated in the bottom panel of Figure [S] (Interestingly, that 
panel also illustrates that, in this case, the mutual information-based measure performed much worse at 
selecting SNPs that could then be used to cluster individuals.) 

At root, statistical leverage provides one way to quantify the notion of "eigenvector localization," i.e., 
the idea that most of the mass of an eigenvector is concentrated on a small number of coordinates. This 
notion arises (often indirectly) in a wide range of scientific computing and data analysis applications; 
and in some of those cases the localization can be meaningfully interpreted in terms of the domain from 
which the data are drawn. To conclude this section, we will briefly discuss these issues more generally, 
with an eye toward forward-looking applications. 

When working with networks or graph-based data, the so-called network values are the eigenvector 
components associated with the largest eigenvalue of the graph adjacency matrix. In many of these 
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K-mGans on all SNPs 





Figure 5: Pictorial illustration of the clustering of individuals from nine populations typed for 9, 160 
SNPs, from [10]. Top panel illustrates fc-means clustering on the full data set. Bottom panel plots, 
as a function of the number of actual SNPs chosen and for three different SNP selection procedures, 
the correlation coefficient between the true and predicted membership of the individuals in the nine 
populations. See the text and fTO for details. 



applications, the network values exhibits very high variability; and thus they have been used in a number of 
contexts |175j , most notably to measure the value or worth of a customer for consumer-based applications 
such as viral marketing [176) . Relatedly, sociologists have used so-called centrality measures to measure 
the importance of individual nodes in a network. The most relevant centrality notions for us are the 
Bonacich centrality 11771 and the random walk centrality of Newman |168| . both of which are variants 
of these network values. In still other network applications, effective resistances (recall the connection 
discussed in Section |6^ have been used to characterize distributed control and estimation problems |178) 
as well as problems of asymptotic space localization in sensor networks |179| . 

In many scientific applications, localized eigenvectors have a very natural interpretation — after all, 
physical particles (as well as photons, phonons, etc.) themselves are localized eigenstates of appropri- 
ate Hamiltonian operators. If the so-called density matrix of a physical system is given by p(r, r ) — 
V'i('')"^V'i('' )? then if V is the matrix whose column vectors are the normalized eigenvectors tpi, i — 
1, . . . , s, for the s occupied states, then P — VV'^ is a projection matrix, and the charge density at a 
point Vi in space is the diagonal element of P. (Here the transpose actually refers to the Hermitian 
conjugate since the wavefunction ipi(r) is a complex quantity.) Thus, the magnitude of this entry as well 
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as other derived quantities like the trace of p give empirically-measurable quantities; see, e.g., Section 6.2 
of |180] . More practically, improved methods for estimating the diagonal of a projection matrix may have 
significant implications for leading to improvements in large-scale numerical computations in scientific 
computing applications such as the density functional theory of many-atom systems |180[ I181j . 

In other physical applications, localized eigenvectors arise when extreme sparsity is coupled with ran- 
domness or quasi-randomness. For example, [1821 1183] describe a model for diffusion in a configuration 
space that combines features of infinite dimensionality and very low connectivity — for readers more fa- 
miliar with the Erdos-Renyi G„p random graph model |184j than with spin glass theory, the parameter 
region of interest in these applications corresponds to the extremely sparse region ^/n <p < log n/n. For 
ensembles of very sparse random matrices, there is a localization-delocalization transition which has been 
studied in detail [1851 11861 1187[ [TS5] . In these applications, as a rule of thumb, eigenvector localization 
occurs when there is some sort of "structural heterogeneity," e.g., the degree (or coordination number) 
of a node is significantly higher or lower than average. 

Many large complex networks that have been popular in recent years, e.g., social and information net- 
works, networks constructed from biological data, networks constructed from financial transactions, etc., 
exhibit similar qualitative properties, largely since these networks are often very sparse and relatively- 
unstructured at large size scales. See, e.g., [1891 11901 1191[ 1192) for detailed discussions and empirical 
evaluations. Depending on whether one is considering the adjacency matrix or the Laplacian matrix, lo- 
calized eigenvectors can correspond to structural inhomogeneities such as very high degree nodes or very 
small cluster-like sets of nodes. In addition, localization is often preserved or modified in characteristic 
ways when a graph is generated by modifying an existing graph in a structured manner; and thus it 
has been used as a diagnostic in certain network applications [1931 1194] . The implications of the algo- 
rithms described in this review remain to be explored for these and other applications where eigenvector 
localization is a significant phenomenon. 

7 A few general thoughts, and a few lessons learned 

7.1 Thoughts about statistical leverage in MMDS applications 

One high-level question raised by the theoretical and empirical results reviewed here is: 

• Why are the statistical leverage scores so nonuniform in many large-scale data analysis applications? 

The answer to this seems to be that, intuitively, in many MMDS application areas, statistical models 
are implicitly assumed based on computational and not statistical considerations. That is, when com- 
putational decisions are made, often with little regard for the statistical properties of the data, they 
carry with them statistical consequences, in the sense that the computations are the "right thing" or the 
"wrong thing" to do for different classes of data. Thus, in these cases, it is not surprising that some 
interesting data points "stick out" relative to obviously inappropriate models. This suggests the use of 
these importance sampling scores as cheap signatures of the "inappropriateness" of a statistical model 
(chosen for algorithmic and not statistical reasons) in large-scale exploratory or diagnostic applications. 
A second high-level question raised by the results reviewed here is: 

• Why should statistical leverage, a traditional concept from regression diagnostics, be useful to 
obtain improved worst-case approximation algorithms for traditional NLA matrix problems? 

Here, the answer seems to be that, intuitively, if a data point has a high leverage score and is not an 
error then it might be a particularly important or informative data point. Since worst-case analysis 
takes the input matrix as given, each row/column is assumed to be reliable, and so worst-case guarantees 
are obtained by focusing effort on the most informative data points. It would be interesting to see if 
this perspective is applicable more generally in the design of matrix and graph algorithms for other 
MMDS applications. 
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7.2 Lessons learned about transferring theory to practice 

More generally, it should be emphasized that the randomized matrix algorithms reviewed here represent 
a real success story in bringing novel theoretical ideas to the solution of practical problems. This often- 
cited (but less frequently-achieved) goal arises in many MMDS applications — in fact, bridging the gap 
between TCS, NLA, and data applications was at the origin of the MMDS meetings |117|ll95(ll96j . which 
address algorithmic and statistical aspects of large-scale data analysis more generally. Not surprisingly, 
the widespread interest in these algorithms is not simply due to the strong worst-case bounds they achieve, 
but also to their usefulness in downstream data analysis and scientific computing applications. Thus, it is 
worth highlighting some of the "lessons learned" that are applicable more generally than to the particular 
algorithms and applications reviewed here. 

• Objective functions versus certificates. TCS is typically concerned with providing bounds on 
objective functions in approximate optimization problems (as in, e.g., and (|16p ) and makes 
no statement about how close the certificate (i.e., the vector or graph achieving that approximate 
solution) is to a/the exact solution of the optimization problem (as in, e.g., ()12|) ). In machine 
learning and data analysis, on the other hand, one is often interested in statements about the quality 
of the certificate, largely since the certificate is often used more generally for other downstream 
applications like clustering or classification. 

• Identifying structure versus washing out structure. TCS is often not interested in identifying 
structure per se, but instead only in exploiting that structure to provide fast algorithms. Thus, 
important structural statements are often buried deep in the analysis of the algorithm. Making 
such structural statements explicit has several benefits: one can obtain improved bounds if the tools 
are more powerful than originally realized (as when relative-error projection algorithms followed 
additive-error projection algorithms and relative-error sampling algorithms simply by performing 
a more sophisticated analysis) ; structural properties can be of independent interest in downstream 
data applications; and it can make it easer to couple to more traditional numerical methods. 

• Side efTects of computational decisions. There are often side effects of computational decisions 
that are at least as important for the success of novel methods as is the original nominal reason for 
the introduction of the new methods. For example, randomness was originally used as a resource 
inside the algorithm to speed up the running time of algorithms on worst-case input. On the other 
hand, using randomness inside the algorithm often leads to improved condition number properties, 
better parallelism properties on modern computational architectures, and better implicit regular- 
ization properties, in which case the approximate answer can be even better than the exact answer 
for downstream applications. 

• Significance of cultural issues. TCS would say that if a randomized algorithm succeeds with con- 
stant probability, say with probability at least 90%, then it can be boosted to hold with probability 
at least 1 — 5, where the dependence on 6 scales as 0{\og{l/S)), using standard methods 128 . Some 
areas would simply say that such an algorithm succeeds with "overwhelming probability" or fails 
with "negligible probability." Still other areas like NLA and scientific computing are more willing 
to embrace randomness if the constants are folded into the algorithm such that the algorithm fails 
with probability less than, say, 10^^^. Perhaps surprisingly, getting beyond such seemingly-minor 
cultural differences has been the main bottleneck to technology transfer such as that reviewed here. 

• Coupling with domain experience. Since new methods almost always perform more poorly 
than well-established methods on traditional metrics, a lot can be gained by coupling with domain 
expertise and traditional machinery. For example, by coupling with traditional iterative methods, 
minor variants of the original randomized algorithms for the LS problem can have their e dependence 
improved from roughly 0(l/e) to 0(log(l/e)). Similarly, since factors of 2 matter for geneticists, by 
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using the leverage scores as a ranking function rather than as an importance samphng distribution, 
greedily keeping, say, 100 SNPs and then filtering to 50 according to a genetic criterion, one often 
does very well in those applications. 

8 Conclusion 

Randomization has had a long history in scientific applications |197[ I198j . For example, originally devel- 
oped to evaluate phase space integrals in liquid-state statistical mechanics, Markov chain Monte Carlo 
techniques are now widely-used in applications as diverse as option valuation in finance, drug design in 
computational chemistry, and Bayesian inference in statistics. Similarly, originally developed to describe 
the energy levels of systems arising in nuclear physics, random matrix theory has found applications in 
areas as diverse as signal processing, finance, multivariate statistics, and number theory. Randomized 
methods have been popular in these and other scientific applications for several reasons: the weakness of 
the assumptions underlying the method permits its broad applicability; the simplicity of these assump- 
tions has permitted a rich body of theoretical work that has fruitfully fed back into applications; due to 
the intuitive connection between the method and hypothesized noise properties in the data; and since 
randomization permits the approximate solution of otherwise impossible-to-solve problems. 

Within the last few decades, randomization has also proven to be useful in a very different way — as a 
powerful resource in TCS for establishing worst-case bounds for a wide range of computational problems. 
That is, in the same way that space and time are valuable resources available to be used judiciously 
by algorithms, it has been discovered that exploiting randomness as an algorithmic resource inside the 
algorithm can lead to better algorithms. Here, "better" typically means faster in worst-case theory when 
compared, e.g., to deterministic algorithms for the same problem; but it could also mean simpler — 
which is of typically interest since simpler algorithms tend to be more amenable to worst-case theoretical 
analysis. Applications of this paradigm have included algorithms for number theoretic problems such as 
primality testing, algorithms for data structure problems such as sorting and order statistics, as well as 
algorithms for a wide range of optimization and graph theoretic problems such as linear programming, 
minimum spanning trees, shortest paths, and minimum cuts. 

Perhaps since its original promise was oversold, and perhaps due to the greater-than-expected difficulty 
in developing high-quality numerically-stable software for scientific computing applications, randomiza- 
tion inside the algorithm for common matrix problems was mostly "banished" from scientific computing 
and NLA in the 1950s. Thus, it is refreshing that within just the last few years, novel algorithmic 
perspectives from TCS have worked their way back to NLA, scientific computing, and scientific data 
analysis. These developments have been driven by large-scale data analysis applications, which place 
very different demands on matrices than traditional scientific computing applications. As with other 
applications of randomization, though, the ideas underlying these developments are simple, powerful, 
and broadly-applicable. 

Several obvious future directions seem particularly promising application areas for this randomized 
matrix algorithm paradigm. 

• Other traditional NLA problems and large-scale optimization. Although least squares 
approximation and low-rank matrix approximation arc fundamental problems that underlie a wide 
range of problems, there are many other problems of interest in NLA — computing polar decomposi- 
tions, eigenvalue decompositions, Cholesky decompositions, etc. In addition, large-scale numerical 
optimization code often uses these primitives many times during the course of a single computa- 
tion. Thus, for example, some of the fast numerical implementations for very overdetermined least 
squares problems that were described in Section [4.51 can in principle be used to accelerate interior- 
point methods for convex optimization and linear programming. Working through the practice in 
realistic computational settings remains an ongoing challenge. 
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• Parallel and distributed computational environments. In many applications, communica- 
tion is more expensive than computation. This is true both for computations involving a single 
machine — recall recent developments in multicore computing — as well as for computations run 

across many machines such as in large distributed data centers. In some cases, such as with 
Gaussian-based random projections, computations can be easily parallelized; and numerical im- 
plementations for both the least squares approximation problem and the low-rank approximation 
problem have already exploited this. Taking advantage of modern computer architectures and 

systems requirements more generally is a substantial challenge. 

• Sparse graphs, sparse matrices, and sparse projections. Sparsity is a ubiquitous property 
of data, and one which is a challenge since vanilla applications of randomized algorithms tend to 

densify the input data. In some cases, sparsity in the input is structiucd and can be exploited by the 
randomized algorithm; while in other cases it is less structured but it can be respected with existing 
projection methods. More generally, sparse projection matrices are of interest — such projection 
matrices make it easy to perform incremental updates in data streaming environments, they can 
make it easier to perform matrix- vector products quickly, etc. Problems of recovering sparse signals 
have been approached by researchers in theoretical computer science, applied mathematics, and 
digital signal processing; and in many cases the approaches are somewhat orthogonal to that of the 
work reviewed here. 

• Laplacian matrices and large informatics graphs. Laplacians are fundamental matrices as- 
sociated with a graph, and they permit many of the randomized matrix algorithms we have been 
discussing to be applied to graph-based data. In some cases, the goal might be to sparsify an input 
graph; but more typically graphs arising in data applications are sparse and irregular. In the case of 
large social and information networks, for example, it is known that, while often there exists good 
small-scale clustering structure, there typically does not exist good large-scale clustering structure. 
Part of this has to do with the heavy-tailed properties in these graphs, which imply that although 
there may exist a small number of most important nodes in the graph, these nodes do not cap- 
ture most of the information in the data. This presents numerous fundamental challenges for the 
algorithms reviewed here, and these challenges have just begun to be addressed. 

• Randomized algorithms and implicit regularization. In many cases, randomized algorithms 
or their output are more robust than their deterministic variants. For example, algorithms may 
empirically be less sensitive to pivot rule decisions; and their output may empirically be "nicer" and 
more "regular" — in the sense of statistical regularization. Existing theory (reviewed here) makes 
precise a sense in which randomized matrix algorithms are not much worse than the corresponding 
deterministic algorithms; but quantifying a sense in which the output of randomized matrix algo- 
rithms is even "better" than the output of the corresponding deterministic algorithms is clearly of 
interest if one is interested in very large-scale applications. 

In closing, it seems worth reminding the reader that researchers often look with great expectation toward 
randomness, as if a nai've application of randomness will somehow solve all of one's problems. By now, 
it should be clear that such hopes are rarely realized in practice. It should also be clear, though, that a 
careful application of randomness — especially when it is coupled closely with domain expertise — provides 
a powerful framework to address a range of matrix-based problems in modern massive data set analysis. 
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